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ABSTRACT 


This report considers the dynamics of spin stabilized 
spacecraft with movable appendages and is an extension of the 
research reported in Part I (May 1974 - May 1975) where two basic 
types of appendages were treated: (1) a hinged type of fixed 

length whose orientation with respect to the main part can vary and 
(2) a telescoping type of varying length which could represent 
extensible antennas or a tether connected to the main part of 
the spacecraft. 

In this report, the dynamics and stability of a spin stabilized 
spacecraft with a hinged appendage system are treated analytically 
and numerically. The hinged system consists of a central hub with 
masses attached to (assumed) massless booms of fixed length whose 
orientation relative to the main part can change. The general three 
dimensional deployment dynamics of such a hinged system is considered 
here without any restriction on the location of the hinge points. 

The equations, of motion for the hinged system, with viscous damping 
at both hinge points, are linearized about the nominal equilibrium 
position where the booms are orthogonal to the nominal spin axis 
for the case of two dimensional and three dimensional motion. 

Analytic stability criteria are obtained from the necessary condition 
on the sign of all the coefficients in the system characteristic equation 
For stability it is found that (hinge) damping must be present and that 
for limiting cases, where the spin axis is an axis of symmetry, certain 
inequalities relating the hinge point offset coordinates to the moment 
of inertia ratio and end masses must be satisfied. 


ii 


The hinge damping is always required for the nominal deployment 
of hinge members. Numerical results indicate that the rate 
damping is required to remove the transverse angular velocities 
effectively since the hinge dampers alone do not provide satis- 
factory nutation time constants for the system parameters selected . 

Next, the control of a spin-stabilized spacecraft with 
movable telescoping appendages is considered with an application 
of the linear regulator problem. It is assumed that the spacecraft 
consists of a rigid central hub and one or two movable telescoping 
booms (with end masses) which are in general linearly offset from 
the hub principal exes. The equations of rotational motion are 
developed and linearized about either of two desired final states: 

(1) a flat spin about only one of the hub principal axes or (2) 
a zero inertial angular velocity state. A control law for the 
boom end mass position is sought such that a quadratic cost functional 
involving the weighted components of angular velocity plus the control 
is minimized when the final time is unspecified. For such a system 
the computation of the control involves the solution of the matrix 
Riccati algebraic equation. For three axis control more than one 
offset boom (orthogonal to each other) is required. When only two-axis 
control is required and a single boom is offset in only one direction, 
an analytic solution of the matrix Riccati equation is achieved; when 
this system is used for reducing the nutation angle of a spinning 
spacecraft the time constant obtained is one order of magnitude smaller 
than previously achieved using non-optimal control logic. For the 
general case of three axis control results are obtained numerically. 


m 


The problem of optimal control with a minimum time criterion 
has been examined analytically for the special case of a single 
offset boom where it is assumed that the initial conditions are 
such that the system can be driven to the equilibrium (rest)- 
state with only a single switching maneuver in the bang-bang 
optimal sequence. For this system it is possible to obtain an 
analytical solution for the switching and final times in terms 
of the initial conditions and magnitude of the maximum value of 
the control force. 
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fit, dimensionless time 
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I. INTRODUCTION 


This report will describe a continuation of the research 
already accomplished during the first year (May 1974-May 1975) 
on the dynamics of spin stabilized spacecraft with movable 
appendages.^ Part I concentrated on the analysis of the motion 
of a spinning spacecraft during the deployment of two types of 
movable appendages - the telescoping rod type of varying length 
during deployment* and fixed length appendages whose orientation 
with respect to the main hub can vary. In addition the use of 
these appendages to detumble a spacecraft with a random spin 
to achieve final states of (1) close to zero inertial angular 
rate and (2) a final spin rate about one of the principal axis 
was also considered.^ In this report (Part II), the following 
topics are treated; the dynamics and an extensive stability 
analysis of a spacecraft with hinged appendages; an examination 
of linear optimal control theory as applied to the deployment 
maneuver of a telescoping boom system by selecting different 
integrand functions in the cost functional; and the time optimal 
control of a nutating spacecraft using a single offset telescoping 
boom system. 

The first phase of the current study will examine the general 
three dimensional motion of a spacecraft with hinged appendages. 

The hinged system consists of a central hub with masses attached 
to (assumed) massless booms of fixed length whose orientation relative 
to the main part can change. 
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The booms are attached at a given radius from the spin axis and 
v/hen the booms are released they swing outward from the spin 
axis under the influence of centrifugal forces. 

The dynamics of this type of fixed length appendage system 

during the deployment maneuver has been previously studied only 

for the case where the transverse components of the angular 

velocity vector are assumed to be zero throughout deployment 

and where the hinge points are located on the hub's principal 
2 

transverse axes. The motion and stability of such a system 
will be studied, analytically for special cases, and numerically 
for the general case. It is assumed that there is no restriction 
on the location of the hinge points. 

The second phase of the study will consider the control of 
a spin-stabilized spacecraft with movable telescoping appendages 
with an application of the linear regulator problem. It is 
assumed that the spacecraft consists of a rigid central hub and 
one or two movable telescoping booms (with -end masses) which are 
in general linearly offset from the hub principal axes. An 
advantage of such a telescoping system as used in the control of 
a tumbling spacecraft is its potential reuse. The booms can be 
retracted at the end of each control sequence (for very fine 
pointing requirements the small residual angular velocity 
components could be removed by temporarily activating on-board 
damping devices). 


Methods of recovering a tumbling spacecraft using three 

sets of telescoping booms deployed along the hub principal 

3 

axes were examined in a recent paper. From an application of 
Lyapunov's second method (using a modified form of the rotational 
kinetic energy as a Lyapunov function) sequences of boom extension 
maneuvers can be determined so that the spacecraft will approach 
either of two desired final states: close to a zero inertial 

angular velocity state, or a final spin rate about only one of 

3 

the principal axes. Although these sequences of boom maneuvers 
will result in recovery of a tumbling system they will not, gen- 
erally, satisfy any criteria of optimal control theory. 

4 5 

In a related problem, Edwards and Kaplan ’ have examined 

the problem of detumbling a spacecraft using a single internal 

mass which is constrained to move along a linear track within the 

vehicle. A control law was selected such that the net average time- 

rate of change of excess rotational kinetic energy will be negative. 

It was concluded that system performance could be improved by making 

the control mass as large as possible and allowing for larger amplitude 

4 5 

motions along the track. ’ The amplitude of any internal device 
would be limited by the size of the spacecraft; on the other hand, 
externally controlled appendages will be subject to external per- 
turbations (such as solar pressure) and if sufficient length is 
extended, the flexibility of such a structure must be considered. 

Recently Amieux and Liegeois^ have applied linear optimal con- 
trol theory to the two-axis control of a spinning spacecraft system 
using a motor controlled internal single-degree-of- freedom pen- 
dulum. 
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It was assumed that the spin rate remains much greater than the 
magnitude of the transverse angular velocity vector during the 
decay of the nutation angle. 

In the present study both tv/o and three axis optimal control 
of a spin-stabilized spacecraft using movable telescoping appendages 
will be considered. The control will be implemented by varying 
the position of boom end masses which are considered large in 
comparison with the boom mass itself. 

The difficulty in determining a control sequence of extension 

rates for different pairs of telescoping booms which would yield 

a time-optimal recovery of a tumbling spacecraft has already been 

reported. The problem has been that when the equations are 

written in standard state form - e.g. for a case of two sets of 

booms parallel to the spin axis - (where symmetry about this axis 

is maintained during extension), the control function (two different 

7 1 

extension rates) is non-linearly coupled with the state variables. ’ 
Here the problem of time optimal control will be examined analytically 
for the special case of a single telescoping offset boom under a 
particular range of initial conditions. 
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II. MOTION AND STABILITY OF THE HINGED SYSTEM 

1 • Derivation of Kinetic Energy 

The hinged system to be studied is shown schematically in 
Fig. 2.1(a). The co-ordinate system representation is shovm in 
Fig. 2.1(b). The system consists of a central hub with masses 
attached to massless booms of constant length which in turn 
are attached to the main spacecraft at radius rp. The end masses 
are released at t = to and thereafter swing out from the spin 
axis. The angles between the booms and the spin axis are denoted 
by and og as shown in Fig. 2.1(b) and are assumed to be zero 
initially. A special case of this type was considered in Ref. 2 
(where it was assumed that the transverse angular velocities during 
deployment remained at zero), but here we consider the general three 
dimensional deployment dynamics. 

The development of the kinetic energy of this type of hinged 
system from first principles is considered below: 

The total kinetic energy of the system, in terms of rotational 
and translational energies, can be written as, 

~ T-t const, due to (circular) orbital motion (2.1) 

where 

Vi + l2«>2 + (2-2) 

1 2 2 
T 4 .= i M + z m. V , 

t 2 M/cm 1 m./cra ( 23 ) 

(M = mass of main body) 
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From the definition of the center of mass of the system: 


cm/o 


m 2 r^-/o 
i=1 ^ 
n 

M + 2 m 


(2.4) 


i=l 


1 


where point 'o' is the center of the coordinate system (and also 
the center of mass of the main hub) and the masses are assumed 

i 

to be equal (m^ = m). The velocity of the various components 
relative to the system center of mass may be expressed: 


\ , = V / + V / 

m^/cm m^-/o o/cm 


(2.5) 


V = V + V 
M/cm M/o o/cm 


( 2 . 6 ) 


The components appearing in Eqs. (2.5) and (2.6) can be further 
represented as: 


'fm^/o "" ^i 


V =0 

o/M 


'^o/cm ^cm/ 0 '^cm/ o 


m 2 r. 


M + 2m. 

1 


(2.7) 

( 2 . 8 ) 

(2.9) 


After substitution of Eqs. (2.7), (2.8) and (2.9) into Eq. (2.3), 
the translational energy may be expressed as 


T = 2 Iv + 7 £ m- IV I 
t , o/cm ^ ^ ^^i/o 


2 1 ' O/cm' 1 m.j^Q o/cm 


(2.10) 
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After some algebraic manipulations, we obtain, 


T = § . VJ - (E V, . z: V.) 

t ' 1 2M ^ ^ 


( 2 . 11 ) 


where 


V. = r, 


+ wx r 


1 


rot 

M = M + z m^ 

Thus, the total kinetic energy of the system is given by: 


1 9 2m 

T = "o (It CO? + + Ioto«) + "p _Z_ (V,. . V,. ) 


2 ^M“l ^2“2 " ^ i=l ' i 1 

!L ( § V. . E VJ + const. 

2M i=1 ’ 1=1 ’ 


( 2 . 12 ) 


As an example, we consider the case from Ref. 2 where m = m/2, 

a-| = ct 2 = a, = I and to^ = 6. The kinetic energy is then obtained 

as (neglecting orbital motion) 


1 


2 . m p 2 .2 .2 


T = -oIe +^[s- a + e (r + 5, sin a) ] 


2 (M + in) “ “ 


(2.13) 


which corresponds identically with Eq. (18) of Ref. 2, which was 
presented without development. 

Next, a more general case of the hinged deployment system con- 
sidered is shown in Fig. 2.2. Here there is no restriction on the 
location of the hinge points. The co-ordinates of the two masses 
are given by 


(2.14) 



= To + sin = “(»'o ^ “ 2 ^ 

2 ^ = a^ - K, cos a^ z^ = a^ - cos a^ 

Here 'a^‘ is the offset of the hinge point from the '2' axis. Upon 

substituting Eqs. (2.14) into Eq. (2.12), and after algebraic simplifi- 
cations* the resulting equation for kinetic energy is. 

1222 
T = 2 [1^03^ + 12^2 

2. 

+ |. [{2(ro+a*+/)+2ji {rQ(sina-j+sina2)-a^(cosa.|+cosa2) )} W] 

2 222.2 

+ {2a^-2a*5.(cosa.j+cosa2)+£- (cos a-|+cos a 2 )> W 2 

+ {2rp+2rQji(sina^+sina2)+5' (sin a-j+sin 02 )} 

- ill {a^;(sinai-sina2)-ro(cosa.j-cosa2)} 

- 2£,2(sinaiC0Sa.j-sina2C0Sa2)} 

+{2i,2(a^ - a2) + {a^ (r^ sina^ - a^cosa^) 

”“2^’^0 ^^*^“2 " ^*C0Sa2)}} + £^(a^ + a|)] 

■2pSr 2)l2cos(o,+ci2)-4a,l(cosa,+cosc2)^ 

+ {2a^- ^{cosa^+COSa^) (sina^-sina^)^ u| 

-2a(sina^-sina2) {2a*-il(cosa.| + cosa^)) WgWg 
+2UJl(a^-a2) + ^ cos (cc]+a2)(“|"“2^ 

-2a^(cOSa.ja.| - COSa2 CI 2 )} til'll Jt^{a|+oe2“2otia2COS(ai+a2) }]+C0nst, (2.15) 
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2. Development of Equations of Motion 

The equations of motion in the five variables: oj^, 

O 

a and are developed using the Quasi-Lagrangian formulation 
for m., i = 1,2,3, and the general Lagrangian formulation for 
the variables , a^. The equations of motion for this system, 
neglecting external torques, can be represented as: 


_d 

3T 

3T 

+a)o 

3T 


0 

(2.16) 

dt 

8u)i 

^3 3(1)2 

2 

8(1)3 




d 

dT 

3T 

+ (!)_ 

3T 


0 

(2.17) 

dt 

3(1)2 

^1 3^3 

3 

3(01 


_d 

3T 

3T 

+(i). 

3T 


0 

(2.18) 

dt 

80)3 

““2 3w| 

1 

3(1)2 



d 

3T 

_ 3T 



0 


(2.19) 

It 

3a^ 

3o(l 

3a' 

1 


d 

9T 

. IL 



0 


(2.20) 

dt 

^“2 

3«2 

da, 

2 




where 

T = Total kinetic energy of the system 
F = Rayleigh dissipation function 

The dissipation function which accounts for linear viscous damping, 
assumed to be present about the hinge points, is given by: 

F = ^ c^ (a^ + 612) (2.21 ) 

where 

c = the hinge damping coefficient 
With the the approximation: m^/M « m or (m/M « 1), the equations 
of motion are obtained as follows: 
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-m{2a*-2a^£(ca^+ca2)+ii2(c2a^+c2a2)-2ro-2rQ^(sa^+sa2)}a)2W3 

+2ms, { rQ{ca^a^ + Ca^a 2 ) + 501202 )^ 

2 2 

+m^ {a* (sa^-sa2)-rp(ca^-ca2) - |- (s2ai-s2a2)}(w3 *• ^ 2 ) 

+in{)l{ai-a2)+(ai)^(roCaj+a^Sai ) + a-j (rgSa^-a^Ca-j ) 

l2a)2-(l3-I])“3‘^l+"’^2a^"2a*Ji(ca^+Ca2)+a^(c a^+c 

-ni5.{a^(sa-j-Sa2) - ro(ca-j-Ca2) - |- (s2a-|-s2a2)} a)^ 

+mJl{a*(sa^“Sa2) - ro(ca^-Ca2) “ I (s2a^-s2a2)} t0ico2 

-m{/(Sa^+S^a2) -2(a^ H-2-^) + 25,a^(ca^+Ca2)>W3W-| 

+tn{2a^S,(sc4^a^+Sa2a2) - £^{s2a^a^ + 5202612)} ^2 

+m£{x(ai~a2) -2a^(ca^6t-j-Ca2a2) + 2.{c2a^6i|-c2a2a2) > ^^=0 

2 2 

I3‘^3"(ll'l2^“l“2 f^^2ro^+2r^il{sa-|+sa2) + i^^(s a^+s 02)} W3 

-mjl{a*(sai-sa2) - ro(ca^-Ca2) - |-(s2a^-s2a2)} 0)2 
-m{2(ro+£2) + 2iro{sa^+sc£2) - S-^ (c‘'ai+c^a2) ) o)-jO)2 
-mjl{a^{sa^-sa 2 ) - rQ(ca^-Ca 2 )- •|(s 2 a^“S 2 a 2 )> (o^o)^ 
-m£{2rp{Sa^a|-Sa2a2) - ^ (c2a^a^-c2a2a2) + A (a^-612)} oy^ 
+m{ 2 rQ)l(ca-ja-j+Ca 26 i 2 ) + (s 2 a|a^+s 2 a 2 a 2 )} ‘^ 3 “® 


( 2 . 22 ) 


(2.23) 


(2.24) 
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A + (£+rQSa^-a^ca^)d)^-(roCa^ + |.s2a^)a)2 "(a^Sa^ - |.s2a^)(02 

-(rpCa-i + a*sa-|)u)^+(a^ca^+roSa^ -Jlc2a^ )w2W2+ ^ “^= 0 

2 2 

« Oj -(«. +r|,sa2-a»ca2)i] -(>'o '=“2 I *^^*^“ 0 " f 

-( r ^ ca ^ "(^*^“2 ''■'’ o'“2 ^ “2 ° ® 


(2.25) 


(2.26) 


where 

sa^=sina^ an(J Ca^ = COSa.j 

3, Two Dimensional Motion Analysis 

a. Small ainplitu<Je analysis about equilibrium state 
(i) linearization of the equations of motion 
The equations of motion for the two (dimensional case with no 
offset are obtained by assuming (jj^=a> 2 =a*= 0 * Also the hinge members 
are assumed to move in-phase (a-]=( 22 =a) and the viscous damping 
about the hinge points is assumed to be absent. Eqs. ( 2 . 22 ) - 
(2.26) then reduce to: 

{I3 + 2m(rQ + 2 roH sina + sin^a))w3+2mJi{2rpCosa a 


+ Z sin2a a}u =0 (2.27) 

O 

Z a "(r^COSa + sin2a) O)^ = 0 (2.28) 

These two dimensional motion equations are linearized about the 
nominal equilibrium state: a = tt/2 and W 3 = (J 2 =noniinal spin) 
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The original coordinates can then be related to the variational 
coordinates by: 



e 


~ Q + 5 

where 

£<< 1 » 6 << Q 

The linearized equations which result are: 

= Q + 5(0) - (2.29) 


.. (ro+^) , „ 

e + — £ e = 0 

When hinge damping is present, Eq. 

form 


(2.30) 

(2.30) has the following 


e 



(^ 0 +^) 
•e + 


e 


0 


(ii) stability criteria 

Eq. (2.31) can be written as 

e+d *e+X^e=0 
a 

where 

d^ = c^m and X=Vl+r“p/i' « 


(2.31) 


(2.32) 


For stability of this second order system the following conditions 
must be satisfied: 

(1) > 0 - Positive damping is required 

(2) < 2X c^ < 2mQ /T+r^TT - which implies a lower bound on 
the magnitude of the damping coefficient. 
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(iii) closed form solutions 

The solution of Eq. (2.30) is given by 

e = C0s(xt+'i') (2.33) 

where is the initial displacement from a = ti/2 and y is the 
initial phase angle. The corresponding solution of Eq. (2.31), 
assuming that the stability condition d <Zx is satisfied, is 
given by 

e = A cos (xt+HO (2.34) 

Where A and v are found from the initial conditions. Eq. (2.33) 
describes the oscillatory nature of the hinge members about the 
equilibrium state while Eq. (2.34) indicates the damping of this 
oscillation. 

(iv) numerical results 

The nonlinear equations of motion describing the deployment 

of the hinged system are programmed for numerical integration using 

9 

the Nova 840 computer. The details of the subroutines used , the 
listing of the program and the computer time required are given in 
the section COMPUTER PROGRAMS (at the end of this report). For 
numerical integration Eqs. (2.22) - (2.26) are used with the following 
system parameters (Fig. 2.2 )^: 

= 8.5 slug-ft2; = 10.5 slug-ft^ 

r^ = 1 ft; = 4 ft; = 0 

m = 0.125 slug; = 4.82 rad/sec 
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Figs. 2.3(a) and (b) simulate the dynamic response of the 
two dimensional motion of the hinged system to an initial pertur- 
bation in the hinge angle of 0.1 radians with hinge damping absent 
and then, present, respectively. For the damped case a hinge 
damping coefficient of c^=0*l Ib/ft/sec is selected. These figures 
verify the closed form analytical results obtained in Eqs. (2.33) 
and (2.34) using the small angle analysis about the equilibrium 
state (a^ = a£ = a = 90°, - 4.82 rad/sec). 


b. Large amplitude analysis 
(i) closed form solutions 

Eqs. (2.27) and (2.28) are used to simulate the large amplitude 
two dimensional motion when there is no damping about the hinge points. 
(Throughout the two dimensional analysis, the hinged members are 
assumed to move in phase and there is no vertical offset of the 
hinge points.) From Eq. (2.27), the closed form solution relating 
the spin rate to the hinge angle, with the initial conditions ^^{O) = 

^2, ot=0, is 


10 . 


(t) = 


l 2 + 2 m(r^+£- sina)^ 


(2.35) 


Here it is observed that ^^(t) attains a maximum value when a=0, n, 
etc., and a minimum when a=Tr/2, 3fT/2,etc. Eq. (2.35) may be substituted 
into Eq. (2.28). This can then be integrated once with respect to time 
under the following assumed initial conditions: a(0) = ot(0) = 0. The 

resulting expression gives the phase plane relationship of a with a 


as, 
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a, ~ 


a. 


(2.36) 



2 

ro+^. sina) i sina(l 2 + 2 mro) 
+ 2m(ro + y. sina)^ 


(ii) numerical results 

The deployment of the system from the position where the hinged 

members are initially parallel to the spin axis (a=0) is simulated 

in Fig. 2.4(a) without damping* and in Fig. 2.4(b) with hinge damping. 

2 

The 'x represents the maximum time simulated by Lang and Honeycutt. 

It is seen that without damping the hinged members exhibit a flapping- 
type motion as momentum is exchanged between the hinge and spin motions. 
The hinge motions are not exactly sinusoidal as shown. Fig. 2.4(a) 
verifies the closed form analytic solution obtained relating spin 
rate with hinge angle. This figure also indicates the maximum and 
minimum spin rate obtained during the flapping motion which verifies 
the closed form solution obtained with the large angle analysis. 

Fig. 2.4(b) shows that, with the hinge damping coefficient selected, 
the system can be fully deployed in about 10 seconds. 

Thus for the case of the two dimensional analysis, closed form 
analytical solutions are obtained for small amplitude oscillations 
about the equilibrium state. Also, from the large amplitude analysis, 
a closed form analytical solution relating the spin rate to the hinge 
angle and the hinge angle rate to the hinge angle are obtained. 

4. Three Dimensional Motion Analysis 

a. Small amplitude analysis about equilibrium state 
(i) linearization of equations of motion 
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The three dimensional equations of motion, Eqs. (2.22) - (2.26), 
are linearized about the nominal equilibrium state: to-j = W2 = 

a-| = og = it/2, 6J3 = The original coordinates are related to 
the variational coordinates by: 

a-j = nr/2 + e-j 

= 11/2 + £2 

lOg + 5 

Based on the assumptions that u-j/fl, 032/fi» and are 

small compared to 1, the equations of motion can be approximated by 

a linear set in which 103 = ft + 6(0) =: ft is a constant. With t = ftt 

as the independent variable, a = w-j/ft and 3 = w2/ft as dependent 
variables, and 

h = -f a^ , g = -j (ro+^)2, k = — (r^+i:) 

2p = cymft, f = 1 +(ro/J?')5 n = (I3/I) - V 

(Here the relationship g = 2kf is to be noted) as nondimensional 

constants, the linear equations become: 


(1+g+h) a' + (n+g-h) 3 + 

k( &|"+e-| ) - k( ® 

(2.37) 

(1+h) 3' - (n-h) a =0 


(2.38) 

e-i" + 2 p ej + f e-] + f a' 

' + f 3 =0 

(2.39) 

+ 2pe^+f 

' - f 3 = 0 

(2.40) 


where primes denote derivatives with respect to t. 
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(ii) stability criteria 

The necessary and sufficient conditions for stability of 
the system are obtained by applying the Routh-Hurwitz criterion 
similar to the procedure of Ref. 10. The stability conditions 
for the system are obtained by setting the characteristic 
determinant, corresponding to Eqs. (2.37) - (2.40), given by 


(l+g+h)s 

(n+g-h) 

k(s2+l) 

-k(s2+l) 

-(n-h) 

(l+h)s 

0 

0 

fs 

f 

s^+2ps+f 

0 

-fs 

-f 

0 

s^+2 ps+f 


equal to zero. Here 's' denotes the characteristic exponent. 

The characteristic determinant is expanded to obtain the 
factored sixth order algebraic equation 

(s^+2 Ps+f)[{(l+g+h) (l+h)s^+(n+g-h)(n-h) }(s^+2PS+f) 

-g(s2+l){(l+h)s^+(n-h)}] = 0 (2.41) 

It is seen that the system characteristic equation separates into 
two factors - a second order factor describing a mode where both 
hinged members move in phase as a unit, and a second factor repre- 
sented by a more complex fourth order polynomial. (Such a separa- 
tion of the system characteristic equation was also observed by 
Auelmann and Lane^^ in studing the stability of a spinning space- 
craft with a ball-in-tube nutation damping system.) From the 
quadratic factor the stability condition is found to be: c^ <2mfiv'THr^7^. 

This condition was also obtained earlier in the analysis of the two 
dimensional hinged system. 
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As the stability, also, depends on the fourth order factor 
in Eq. (2.41), we will consider different configurations of the 
system in the following special cases. Case 1(a): 1 3 maximum 
moment of inertia (n> 0 ) with no offset of hinge points (h= 0 ) . 

The fourth order factor in Eq. (2.41) with h=0 reduces to 

p^s*+ + Pis3’+ P2S^ + p^s + p^ = 0 (2.42) 

where 

= 1, p^ = 2P(l+g) 

?2 = ( 1 +g) f -g(l+n) 

P 3 = 2 P(n+g) n, p^ = (n+g) nf -gn 
The non-trivial Routh-Humitz stability conditions are 

P1P2 - ° 

(P]P2 - P0P3) P3 “ P1P4 > 0 ( 2 . 44 ) 

Expansion of Ineqs. (2.43)and (2.44) results in complex algebraic 
relationships involving the system parameters. An attempt to 
algebraically rearrange the terms in Ineqs. (2.43) and (2.44) did 
not yield simplified results. However, it can be seen from con- 
sideration of the signs of each of the coefficients in Eq. (2.42) 
that P|>0 and P 2 > 0 . Since both g and n are positive, therefore, for 
stability, p>0, implying the necessity of (positive) hinge damping. 
Case 1(b); I 3 maximum moment of inertia (n>0) with offset of hinge 
points (hi^O) . 
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The fourth order factor in Eq. (2.41) vnth h f 0 can be 
written as 

+ (qif-q-j) = 0 (2.45) 

I O 

where 

= ( 1 +g+h) ( 1 +h) 
q^ = (n+g-h) (n-h) 

= g(l+h) 
q 0 = g(n-h) 

The necessary and sufficient conditions are obtained from the rule 
of signs of the coefficients in the characteristic equation and also 
the Routh-Hurwitz criterion and can be expressed as: 


la 1 < / 3 

(2.46) 

1 * 1 - V 2 ni 

W^W2 - WqW3 > 0 

(2.47) 

(WqWg - WqW 3 ) W 3 - >0 

(2.48) 


where 

Wo = qo - qz 
Wir 2«o 

W2 = fqo + qi - q2 - qj 
VI 3 = 2 pq^ 

W4 = qif - qj 


- 19 - 


In Ineq. (2.46) the magnitude of the hinge offset from the (1,2) 
hub plane is limited by the differences in the hub moments of 
inertia and the size of the end masses. Case 2(a): 1 3 minimum 

moment of inertia (n<0) with no offset of hinge points (h=0) . 

The necessary condition from the rule of signs of the 
coefficients is obtained as 


(ro+s<) < minimum of 


^(I-Ij) /2m' 


(I-I3) /2mro 


(2.49) 


Case 2(b): I3 minimum moment of inertia (n<0) with offset of hinge 

points (h40) . 

Similar to case 2(a), the necessary condition is obtained 
as 


(>'o'^^) minimum of 


/ 


( I-I3 ) 

2m 


+ a* 


(2.50) 


( I-I3 ) + ^ 


2mr 


It is seen that when a^=0, Ineqs. (2.50) reduce to Ineqs. (2.49). 

The necessary and sufficient conditions for the cases 2(a) and 2(b) 
are obtained in a similar procedure by making n<0 in the corresponding 
inequalities, Ineqs, (2.43), (2.44), (2.47), and (2.48). 


(iii) numerical results 

An example of the three dimensional hinged system dynamics is 
simulated in Figs. 2.5 (undamped) and 2.6 (with hinge damping). 
Initial perturbations in both hinge angles and one of the transverse 
angular velocities are assumed. 
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The oscillatory nature of the system motion about the equilibrium 
state when there is no hinge damping present is seen in Figs. 2.5(a) 
and 2.5(b). Although the hinge damping is effective in reducing 
the amplitudes of the hinge motion (Fig. 2. 6 (a)), the time constants 
associated with the nutation angle decay (Fig. 2.6(b)) are extremely 
long. It is clear that an additional form of nutation damping must 

I 

be added for effective removal of excessive transverse rates. 

The effect of rate damping about the transverse axes to obtain 
favorable nutation decay time constants is shown in Fig. 2.7. The 
rate damping torques about the and ' 2 ' axes are assumed to be 
R^u-i and R^W 2 » respectively. The response of the transverse angular 
rates for rate damping coefficients of Rjj= 1.0 Ib-ft-sec and 2.0 
Ib-ft-sec, respectively, are considered. The simulation results 
(Fig. 2.7) indicate that for effective removal of excessive transverse 
rates rapidly the magnitude of the rate damping coefficient must be 
large. The simulation study shows that the rate damping does not 
have any effect on the hinge angles and spin rate responses which 
are essentially unchanged from Fig. 2.6 for the parameters considered 
here. 

The deployment dynamics of the system when the vertical hinge 
points are offset (a^j^O) from the '2' axis is considered next. The 
stability condition for a* when I 3 is a maximum moment of inertia 
(n>0) is obtained from Ineq. (2.46). With the system parameters 
selected, the limit of offset from Ineq. (2.46) is a^ < ft. 


- 21 - 



The time response of the transverse angular rates with a^ = 1.0 ft 
(< ft) is compared with the deployment dynamics when there is 
no vertical offset (a*=0) in Fig. 2.8. Here it is found that 
the time period of the damped response of the transverse angular 
rates for the case with a vertical offset is larger than the 
period of this response when there is no offset. This may be 
attributed to the redistribution of moments'of inertia due to 
the offset. Also, it is to be noted that the vertical offset 
(a^=1.0 ft) does not have any effect on the time responses of the 
hinge angles and spin rate (essentially the same as shown in Fig. 2.6). 
The magnitude of the rate damping coefficient is assumed to be R^=2.0 
Tb-ft-sec. 

From the small amplitude analysis of the three dimensional 
motion of the hinged system we conclude the- fol lowing: 

1. Hinge damping is required for the nominal deployment of 
hinge members. 

2. For stability, from the rule of signs, certain inequalities 
relating the hinge point offset .to the ratio of the moments 
of inertia and end masses must be satisfied. 

3. Since the hinge dampers alone do not provide satisfactory 
nutation time constants other types of dampers which provide 
direct damping of the transverse angular rates are required, 
and the magnitude of the rate damping coefficient must be 
large for rapid removal of the transverse rates. 
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4. The vertical offset of the hinge points increases 

the time period of the response of the transverse rates, 
which may be due to the redistribution of the moiuents 
of inertia. 

5. The simulation of the deployment dynamics shows that the 
rate damping and vertical offset do not affect the responses 
of the hinge angles and the spin rate with the system para- 
meters and the initial conditions selected here. 

h. Large amplitude analysis 
(i) numerical results 

The general three dimensional motion analysis for large amplitude 
(a-j and are physically free to vary between 0 and 180° - see 
Fig. 2.1(a))is considered in this section. The nonlinear equations 
given earlier (Eqs. (2.22)-(2.26) ) are used for numerical simulation. 

The results are illustrated in Figs. 2.9(a)-(c) for different para- 
meters. From the simulation of the nonlinear system motion, (Figs, 2.9 
(a)-(c)), the observations made are: 

Case (1) : For the values a*=0, c^=0 and R(j=0, the hinge members 

would intersect the hub structure (|al > 180°) as indicated. (It should 
also be noted that with these parameters, the motion of the system 
when linearized about a^=ir/2, w^=0, i^l,2, would be unstable in 
the sense of Routh-Hunvitz) . 

Case (2) : For the values a^=0, =0.1 Ib/ft/sec and R^=2 lb/ ft/sec, 

the hinge motion lies within the physical boundary as indicated. 
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This Case(2), shows with the hinge and rate damping coefficients 
selected, that the system can be fully deployed in about 10 secs 
(Fig. 2.9(a)). The angular rates about the transverse axes are 
nearly removed (Fig. 2.9(b)) and the spin rate reaches a steady- 
state value of 4.1 rad/sec from an initial value of 0^3(0) = 4.82 rad/sec. 
(Fig. 2.9(c)). 

Case (3) : For the values: a*=1.0 ft, c^=0.1 and R(j=2, the 

hinge motion would interfere with the main satellite structure, as 
in Case (1). 

Case (4) : Since there are no criteria to determine the magnitude 

of offset of the hinge points for general nonlinear motion, a value 
of a^=0,5 ft, which is less than 1.0 ft used in Case (3), is selected 
keeping the same values for c^ and as in Case (3). The response 
of the system for these parameters is indicated. The behavior here 
is similar to that shown in Case (2). 

The nutation angle e is defined by the equation 

I /I? + COo 

tane = -7-^ ^ 

'3 “3 

The time response of the nutation angle is shown in Fig. 2.9(c). The 
cases of a^=0 and a*=0.5 ft are considered, The nutation angle reaches 
a maximum of 18° during the deployment and declines to a value of 
less than 1° within 10 secs. 
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The large amplitude analysis of the three dimensional motion 
for the general case of deployment reveals the following: 

1. Hinge damping must always be present for the nominal 
deployment of the hinge members. 

2. Rate damping is needed for the effective removal of the 
transverse angular rates. 

3. The selection of the offset of the hinge points can not 
be done easily as there is no criteria which assures that 
the hinge motion will remain withia physical limits 
(i.e. (aj < it). 


- 25 - 





3 



Xl = 0 X2 = 0 

yi = ^0 + ^ sin «i y2 = -(ro + I sin ^z) 

2i = - £ cos Z2 = t cos «2 


FIG. 2.2. MORE GENERAL CASE OF HINGED DEPLOYMENT SYSTEM. 
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III. EQUATIONS OF MOTION WITH 
TELESCOPIC TYPE CONTROL BOOMS 

The complete equations of motion of a rigid spacecraft with 
telescopic type control booms are developed. The system is assumed 

to consist of (Fig. 3.1) a central hub v/ith center of mass at point 
Q and one or two extendible telescoping booms with end masses m-j and 
m 2 » respectively. The mass along the boom lengths is assumed negli- 
gible in comparison with the end masses. It is assumed that, in 
general, the two booms will be offset from the hub principal axes 
with the coordinates a,b,c,d indicating the amount of offset. 

The generalized vector equation of motion for such a system 

12 

containing a central hub and moving connected masses can be written : 

_ j- N _ 

'■'Q “ 4 "'k (>'k/q) '■q (3.1) 

where Hq refers to the external moments, Q refers to the reference 

• • 

point which is assumed to be at the center of mass of the hub, Fq 
is the inertial acceleration of the reference point and is the 
position vector of mass, m[^, with respect to point Q (Fig. 3.2). 

It should also be noted that is the positi on .vector of the composite 
system center of mass whose position v/ill change with the movement of 
my and The composite c.m. is assumed to move in a circular orbit, 
and it is assumed that coupling between orbital (translational) motion 
and the attitude dynamics is a higher order effect. 


- 41 - 


The angular momentum of the system measured with respect to 
point, Q, has three components, 

" ^b/Q ^m-|/Q ■*' ^m2/Q (3.2) 

where L|^^q describes the momentum of the hub, and Lm^/q describes the 

momentum of mass m^. The hub momentum may be expressed in terms of 

the hub principal moments of inertia and angular velocity components 

as : _ _ 

Lb/q = ^ ^3^3*^ (3.3) 

where T,J,F are unit vectors along the hub principal axes, and 

4 

^/Q = ^ i i = (3.4) 

where 7 ^ describes the position of m^. relative to Q(»'^/q). 

We will now consider the inertial acceleration of the reference 
point (Fig. 3.2) 


••q = Rc - '•c 


(3.5) 


In the absence of external perturbation forces and within the previous 


assumptions R = 0, and. 


'•q “ - ''c 


(3.6) 


From the definition of the system center of mass we can relate 


• 9 

7r 




^ M + Em 

where M represents the hub mass and 


Em = m.j + m^ 


(3.7) 
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After substituting Eqs. (3.3), (3.4) and (3.7) into Eq. (3.1) 
and assuming the external torques vanish, the following vector 

rotational equation results: 

♦ * • * * 

+ Us {F-| X + ^2 ^ ^1^ 


where 

= m^ (M + ni2)/(M + Em) 

P2 = ni2 (M + m^ )/(M + Em) 

P2 = - m^m2/(M + Em) 

Eq. (3.8) is then expanded using the familiar relationship. 


— — T- 

dt body 

(3.9) 

and for the specific geometry of Fig. 3.1, 


r-j = aT + bj + zF 

(3.10) 

r 2 = xT + cj + dF 

(3.11) 


The acceleration terms r^* (i = 1,2) may be calculated by using^^, 

• • 

j' . “ CO X ( (0 X ) "^ CO x r . 2co x C r^ 1 

• • 

l^*^i^body 


together with Eqs. (3.10) and (3.11). 

The complete nonlinear equations of motion are obtained by expansion 
of Eq. (3.8) and are expressed as 

^l“l CO2CO2 + u-j [(b^+z^) oil “ abco 2 

-azco3 “ ^^^1^2 “2^3 ^ 

+2zz CO1 + bz (to^^ “ ^2^^ ^ + VI2 ^ ^2. 
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+cxojio)2 - 2cx(i)2 " Zdxoj-j + cd 

+P3[2(bc+d2) - (ac+bx) ^2 - (ad+xz) ^2 

-(ad+zx)aj-|W2 + 2(bc-dz) 032^^3 (^c+bx) w-iw^ 

+2dztoi - 2bxo)2 ~ ^xzw^ + (bd+cz) (oj 2^-W2^) + cz] = 0 

• • • 

12^2 + ^3^1 ■ 1^] [abco-j - (a^+z^) W2 

+bzo)3-bzuiW2 atw^w^ + ‘(a^-z^) 

-2zzw 2 az (003^-0)1 2) + az] - y2fcxoj^ 

-(d^+x^)t02 ^*^^3 " cdu)ia32 + CXW2W2 

+ (x2-d^)o32Wi - 2 xxo) 2 + dx(oj2^-WT^) - dx] 

• • * 

~y3C (ac-i-bx)wi -2(ax+dz)w2 + (bd+cz)oj2 

-{bd+cz)to-|W2 + (ac+bx)oj2f^3 + 2(ax-dz)oj20Ji 

-2{ax+dz)w2 + (ad+xz) (0)3^-031 2) 


+XZ-ZX] = 0 


I2W2 + “1^2 " ■*■ 


-(a^+b^) 032 ^ (b^-a^} 03-jO32 
+ 2 azo 3 -j + 2 bzo 32 + ab(a 3 ^ 2-032 


- + bzoo^o), 

2)3 ^ 


+cdo32 - (c^+x^) 032 + (c2-x2) 03-j0)2 

t 

“dxa)2W2 + cdo3203^ - 2xxo32 + cx(w-|2 
-032'^) + cx] - P3 [(ad+xz) w-| + (bd+cz) 

-2(be+ax) 033+2 (bc-ax) 03^ o)^ - (ad+xz) 03203^ 

* • * 

+( bd+cz) 0)3031 + 2 xzo3i ~ ^^^“3 

+ (ac+bx) (o)i^ - 0)2^) + bx] = 0 


(3.13) 


(3.14) 


(3.15) 



Eqs. (3.13) - (3.15) are three coupled, non-linear differential 
equations for the spacecraft dynamics in terms of the angular rates 

, the positions of the control boom end - masses (a,b,z) 

« • 

and (x,c,d), and the corresponding velocities (z,x), and accelerations 
(z,x). These equations are valid irrespective of the physical mechanism 
which causes the control booms to excute their motions. For the 
special case v/here m 2 = 0, these equations correspond identically with 
Eqs, (2) - (4) of Ref. when Ir-] has been defined according to 
Eq. (3. 10). 
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IV. OPTIMAL CONTROL WITH QUADRATIC 
PERFORMANCE CRITERIA 


1. Linearization of Equations of Motion 

The equations of motion developed in Chapter III represented 
by Eqs. (3,13)“(3,15) are linearized about the follovnng desired 
final state; w-j = ug = 0» ^2 = n. It is assumed that w^/n<<! 

(i = 1 , 2 ), |w 2 - 1 , and that and 2 / 1 ^ are both « 1 , 

where is the maximum boom length. In order to develop the 
linear system equations in non-dimensional form we further let: 
a = w-j/n; 3 = Wp/o; 1 + Y “ ^ 2 /'^’ ^i “ ~ 1,2,3) 

where y=reference mass; ^ - x/S-^; c = = a/£^; Cp = 

C 3 = c/P-jy,; denote derivatives with respect to the 

nondimensional time, t, by primes ('). The variational coordinates 
are a, 3 , y whereas c and c represent the control variables and 
describe the end mass positions. 

For the case of a single-boom system (Fig. 4.1), m^ = ^3 " 0, 
and letting u = u-j the linear system can be expressed by: 


h " ^2^ 

-CiCp 


I +c ^ 
^2 ^1 

0 

0 


-03-i 


C^jCj 

0 

0 

Vhere 

CU]=[;" 

+ 0 


0 


a' 

0 


3' 




I3+Ci2+C2^ 


y' 



•* 



m - 



:a 



3 



Y 



m « 

L 



(4.1) 

(4.2) 
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It should be noted that for such a system three-axis control 
with a single boom can not be achieved since y' = 0. Eq. (4.1) 
can be reduced to the standard from: X' = AX + BU, where x"^ = 

(a 3 ) and the elements of A are: 


A-j .j = c-j C2( 1 2 “ ^ 1 ” 1 2 ^ 

A21 = {(l3-i-,)(Ii+C22) + IiC.| 2 )/^ 

A 22 = -c^C2(i3-Ii-l2)/A (4.3) 

’I2C2/A 

and A = I^l2 + ^ 2 ^ ^2 




It is further noted that offset of the single m-| boom from at least 
one of the hub principal axes is required for control, since if c-j = 
C 2 = 0 the system is uncontrolled. 

A final state of zero inertial angular velocity can also be 
considered' by nondimensional izing using an arbitrary reference value 


of oj|^ such that «1 for i = 1,2,3 in place of fi. Here the 

non-dimensional variational coordinates are: a = 6 = 

6 = 0 ) 3 / 101 ^; and other parameters are defined as earlier. The resulting 
linear system reduces to 

X' = A 6(0) X + BU (4.5) 
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where 


6(0) = ai^CO) 

U = ?" + 6(0)c 


Eq. (4.5) is similar to the standard form X' = AX + BU. 

For the two-boom case when d = = 0 and the equilibrium 

state is; w-, = o )2 = possible to reduce 

the linear system to the standard form: 


DX' = EX + FU 


(4.6) 


where the non- zero elements of the D and E matrices are: 

Dv) = + (ui/u) ^ ^2^3 


^12 ^ '^^ 1 ^ 1^2 ^ ^ 3 ^ 1 ^ 3 ^ ~ ^21 

^22 "" ^2 

D33 =13+ (ui/y) + c^^) + (^2/^) ^3^ + 2 (u3/u) C2C3 


Ell ~ ■'^^1^1^2 ^ 3 ^ 1 ^ 3 ^ 

Ei2 = -(I3 " ^2 *^2^ (1^2^^^ ^ 3 ^ ^ ^^ 3 ^^^ ^ 2 ^ 3 ^ 

^21 = ^3 - ^1 "1' 

^22 ~ ^1*“2 ^ ^1^2 


and; 


F = 


U = 


“{^1^2 ” ^3^3^ 

(uiCi) /y 

0 


C" + C 
5 " + fi e* 


0 

0 


(n2C3 + V 3 C 2 ) 


(4.7) 
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(4.8) 



v/here 


f, =-Zc^v^Hv^c^ + 

Eq. (4.6) can be written as X' = AX + BU, v/here 

A = D'^E (4.9) 

B = D‘V (4.10) 

» 

The two-boom offset system will now provide three-axis control. The 
systems, Eqs. (4.1) and (4.6), have been linearized about a desired 
final state of spin about the '3' axis; while Eq. (4.5) has been 
linearized about a final state of zero inertial angular velocity. 

2. Application of the Linear Regulator Problem 

It is desired to design a constant gain regulator for the time 
invariant system 

X' = AX + BU (4.11) 

which minimizes the quadratic performance index 

J = /“ (X% + U^RU)dt (4.12) 

where Q and R are constant positive definite symmetric weighting 

13 

matrices. From the theory of optimal control it is known that 
the control law is of the form: 

U = (4.13) 

where K is the symmetric positive definite solution of the algebraic 
matrix Riccati equation: 

KA + A^K - KBR’^B^K + Q = 0 (4.14) 

In general the solution of Eq. (4.14) must be done using numerical 
algorithms because of the algebraic complexity of the problem. 
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For 

the applications 

considered 

here the 

control will 

be 1 

of the form: 






U = 

a + 6 

= 

+ 5 


(4,15) 

for 

the 

two-axis system; 




and 


hii c + h,2 

3 + 

h,3 T = c" . 


(4.16) 


^2 = 

''21 “ ■*' '’22 

3 +' 

h23 Y = 5“ + 


(4.17) 


for the three-axis system, where the g's and h's represent the 
calculated gain constants. For actual implementation of the control 
it will be necessary to measure the components of the main hub angular 
velocity vector (v;hich are proportional to a,3,y)> perhaps with rate 
gyroscopes, and also to determine the length of z boom extension and 
time rate of x-boom extension, perhaps optically or with a counter 
attached to the motor which drives the booms. The necessary accelera- 
tion terms c" and can then be provided by adjusting the motor 
torque to a sufficient level. In this report it will be assumed that 
all measurements occur instantaneously and that there are no errors 
in measurements. It is hoped to treat these important effects in 
a later study with an application of the estimator problem. 

3. Numerical Results 

a. Two-axis control using a single boom to reduce nutation 
angle 

When a single boom is offset from the z axis and the hub is 
symmetrical (I = Ii = I^) it is possible to solve the matrix Riccati 
equation analytically. 
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It is assumed that the boom is offset only along the axis 
(b = 0) without any further loss of generality. The standard 
form of the equations of motion can be written, from Eqs, (4.1) - 
(4.4), as 


'o 1 + U (4.18) 

n 

»• . 

where 

d = (I3 - i + / (i + 



e = (I3 - i) n 


n = c^/ (I + c^^) 

For this application of two-axis control it is logical to select the 
weighting matrix, Q, in the performance index, J, Eq. (4.12), to have 
diagonal elements inversely proportional to the maximum expected 
value of transverse rate 

f 0 where f = ^ (4.19) 

Q = max . 

O f 

■ J 

and we will select R as a unit weighting matrix (following Ref. 6). 

The expansion of the matrix Riccati equation Eq. (4.14), with 
the use of Eqs. (4.18) and 4.19), yields 

2K^2^ - n^K^^ + f = 0 (4.20) 

-K^ie + 1<22^ - n^K^2^22 " ^ 

- 2K ^ 2 ® ” "^^2 ^ ^ 
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Eqs. (4.20) - (4.22) can then be solved for the elements of the 
two dimensional symmetric 1 ( matrix as: 


h2 = 

[d + /W~ 

+ n^ f 

]/n^ 


(4.23) 

*^22 

+ / f - 2 e 

Ki 2 /n 



(4.24) 

hi = 

K 22 (d ~n^ 

K 12 ) /e 



(4.25) 

where the 

choices of 

signs in 

front of 

the radicals 

are selected 

that K is 

positive definite. 

From Eq, 

(4.13) it is 

seen that the 

control has the form 





U = - 

[nK^2 

(X) 



(4.26) 


As an illustrative example the system parameters and initial 
conditions are selected from Ref, 5, pp. 72-79, where the movable 
mass system is considered as a two-axis nutation damper for the NASA 
21 Man Space Station (Fig. 4.2), 

I = 1.42 X 10^ kg ~m^ (10.5 x 10^ slug-vt^) 

I 3 = 2.03 x 10^ kg -m^ (15.0 x 10^ slug-ft^) 

M = 6.21 x 10^ kg (1.37 x 10^ Ibm) 

0 = 0.314 rad/sec (3 rpm) 

ni = 816 kg (1800 Ibm) 

a =19.8 m (65 ft), b = 0, £ = z = 5.4m (17.72 ft) 

in max 

co,(0) = 0,0391 rad/sec, 0 ) 2 ( 0 ) = 0, (o -- 0.04 rad/sec. 

' 'max 
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This results in the following values for the parameters: 
d = 0,666, e = 0.428, n = 0.113, f = 61.63 
and it is seen that K is positive definite if = 219,13, K ^2 " 
-34.69, Kg 2 = 84.57, For ^ 2 ( 0 ) = 3(0) = 0,, the response of the 


state vector components is then calculated as 

a(r) = 0.428t + 1.262 sin 0.428x)a(0) (4.27) 

3(t) = 2.591 sin (0.428t) a(0) (4.28) 

under the optimal control: 

U (t) = 3.92a - 9.573 = C" + C (4.29) 

The solution for ?(x), assuming c(0) = ?'(0) = 0, is: 

c(t) = 0.655 sin t + 0.417 cos t - e"°*^^'^(0,417 cos 0,428t 
+ 2.056 sin 0.428t) (4.30) 


indicating that the boom will experience a steady state oscillation 
after the initial transient. The initial nutation angle is calculated 
to be 5,0 degrees and for small angles may be approximated by: 

9 + 3^ VI 3 (4.31) 

The decay of nutation angle for this case is shown in Fig. 4.3 
(the curve labeled Q-j for = 17.72 ft.) and compared with two of 
the results from Ref, 5. When the maximum amplitude of boom length 
is 5.4m (17.72 ft.) it is seen that the time constant associated 
with the nutation decay is approximately one order of magnitude 
better when the linear quadratic performance criterion is used instead 
of the non-optimal methods of Refs, 4 and 5. In Fig. 4.4(a) for the 
same case the motion of the boom end mass during nutation decay is 
illustrated. 
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It is seen that v/ithin 30-35 secs, the transient part of this 
motion is removed leaving a remaining steady state oscillation 
with an a amplitude of 13.04 ft as predicted by Eq. (4.30). Also, 
the initial disturbances in both the transverse angular velocity 
components are removed by the motion of the control boom. At 
this point the nutation angle reaches a (small) value within 
mission tolerances and the z boom could be withdrawn to the zero- 
length position and ready for re-use. This operation would require 
a temporary boom-motor command or the activation of internal damping 
in the boom extension mechanism so that the steady state amplitude 
would eventually approach zero. 

Also shown in Fig. 4.3 are the affects of varying the Q 
weighting matrix. For the case already discussed the Q matrix has 
the two dimensional form: 

t 61.63 0 “ 

0 61.63 

which yields the optimal control law given in Eq. (4.29) and repeated 
here: 

U (t) = 3.92a -9.573 (4.32) 

If for the same boom offset and we let the maximum expected value 

of transverse angular velocity be, = 0.1 rad/sec, then 

'max 



with the result that 

U (t) = 0.787a -3.9743 (4.33) 
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The nutation angle decay for this case has a slightly longer time 
constant which can possibly be explained by the smaller negative 
coefficient of B in Eq. (4.33). The dynamic response of this case 
is illustrated in Fig. 4.4(b). The amplitude of the steady state 
boom motion is 7.45 ft which is smaller than the amplitude of 
13.04 ft. for the case considered in Fig. 4.4(a). 

For the same offset' and a reduced z = 3m (9.84ft) it is 

max 

seen that with the control law resulting from Q-j the boom end mass 

will actually exceed a displacement amplitude of 3m and the dynamics 

of such a discontinuous system can not be approached analytically. 

For lOj =0.1 rad/sec. Q 2 yields a control law of the form: 
max 

U (t) = 0.4566a -4.00773 (4.34) 

and the corresponding decay of nutation angle is shown in Fig. 4.3. 

The time response of the case considered here is shown in Fig. 4.4(c). 
The steady state oscillation of the control boom end mass has an 
amplitude of 4.40 ft. and the nutation angle becomes zero after 130 
secs. 

In all the cases considered to this point the spin rate is assumed 
to be constant ( 102 ( 0 ) = fi). Next, we consider the case where the 
initial angular velocity about the '3' axis is not the same as the 
final desired spin rate. The control gains are calculated from the 
linearized equations and the dynamic response of the system is obtained 
from numerical integration of the nonlinear equations of motion of the 
system Eqs. (3.13) - (3.15). 
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9 

The details of the numerical integration computer program listings 

are given in the section COMPUTER PROGRAMS. 

The results of this numerical integration are shown in Figs. 

4.5(a) and (b) with weighting matrices selected as Q 2 and Q^, 

respectively. The maximum amplitude of the boom end mass is assumed 

to be 100 ft. The numerical results indicate that more than 200 secs, 

would be required to achieve a zero nutation angle. The transverse 

angular velocity responses are more oscillatory than those shown 

earlier (e.g. Fig, 4.4); and the angular velocity about the '3' 

axis reaches a steady state value equal to the desired final spin 

rate within 150 secs. For this application, the three-axis control 

is achieved using only a single boom due to the nonlinear coupling. 

As Edwards and Kaplan ’ have previously noted, performance is 

improved by increased boom lengths but in all cases considered here 

where the boom motion remains continuous there is a marked improvement 

in system performance by using the quadratic performance criterion 

from linear optimal control theory. 

b. Three-axis control using two offset booms 

For more complicated applications of two-axis control (a f 0, 

b / 0) with asymmetrical hubs and for the general case of three-axis 

control numerical methods were used to solve the matrix Riccati 
14 

equation . Ref. 14 contains both a sample problem and computer 
listings which were adopted for use with the NOVA 840 computer system. 
For the general case, the dynamics of the controlled system were 
simulated by numerical integration of the general (nonlinear) form 
of Eqs, (3.13) - (3.15), as developed in cartesian components. 
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The computer listings and the program details are given in the 
section COMPUTER PROGRAMS. 

An example of the application of a two-boom system for three- 
axis control is shown in Fig. 4.6(a) where initial perturbations 
are assumed in both the transverse angular velocity magnitude as 
well as the component along the spin ('3') axis. Furthermore the 
initial disturbance amplitude is assumed to be less than 0.04 rad/sec 
in both transverse and spin axis components. The non-zero elements 
of the Q matrix are calculated to be 

f = = (0.314)2/(0.1)2 = 9.86 

The hub moments of inertia, boom end masses and hub mass are the 
same as in Sect. a. Based on this information the matrix Riccati 
equation was solved numerically yielding the following control law: 

= 2.2242a -3.84893 = + 4 (4.35) 

U 2 =-3.1401y = 6" + C (4.36) 

The steady state values of the feedback gains were obtained 
by numerically integrating the matrix Riccati differential equation 
until steady state conditions were obtained. This information was 
then incorporated into the simulation of the general rionl inear 
equations of motion, Eqs. (3.13) - (3.15). It is seen from Fig. 
4.6(a) that within 50 seconds the initial disturbance in both the 
transverse angular velocity component as well as (about the 
spin axis) are removed by the motions of the two control booms. 

After this time the z boom continues to exhibit a steady state 
oscillation according to the control law whereas the x boom 
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reaches a constant steady state amplitude of -47,5 ft (14,48m). 

The offset of the z and x booms are assumed to be equal (65 ft.) 

Fig. 4.6(b) illustrates the dynamics of the system for a change 
in the offset value of the x boom from 65 ft. to 45 ft. The 
responses in Fig. 4.6(b) are similar to those shown in Fig. 4.6(a) 
except that the constant steady state amplitude of x is reduced . 
to a value of -43.4 ft. • 

In an actual mission after the disturbances have been removed 
these steady state motions could be stopped by physically retracting 
the booms to the zero position. During all boom maneuvers simulated 
here the extension (or retraction rates) fall within O-lOft/sec which 
is thought to be within that currently obtainable with such mechanisms 

Another example of three-axis control using two offset booms 
is depicted in Fig. 4.7. The initial conditions are identical to 
those of Fig. 4.6(a) as are all parameters relating to maximum boom 
lengths, masses, and offsets. In Fig. 4.7 it is assumed that the 
main part of the spacecraft is no longer symmetric with 

I^ = 1.42 X lO^kg-nf (10.5 x 10^ slug-ft^) 
and 

I^ = 1.69 X 10^kg-m2 (12.5 x 10^ slug-ft^). 

It is seen that for the same weighting of the state variables and 
control, the system in Fig. 4.7 requires a longer settling time to 
reach steady state and that the steady-state amplitude of the x-boom 
is slightly reduced to about -45ft. 
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c. Application of offset system during the terminal phase 


of detuinbling maneuver. 

A recent study has shown that a tumbling spacecraft may 
be recovered by extending three sets of telescoping booms along 
the hub principal axes. With an application of Lyapunov's second 
method^ sequences of boom extension maneuvers can be determined 
so that the spacecraft will approach either of two desired final 
states: close to a zero inertial angular velocity state or a final 

spin about only one of the principal axes. The numerical results 
of Ref. 3 show that there is always a small residual remaining in 
the cG'.iiponents of the state vector and that extremely long times 
(and boom lengths along the principal axes) are required for a 
significant reduction in the residual nutation angle when a final 
spin is the desired final state. It is suggested that the offset 
telescoping boom system can be used as a means of rapidly reducing 
the residual components of the state vector after an initial detumbling 
maneuver using the strategy of Ref. 3, 

As an example of the application of such a system, we consider 


•the NASA 21 Man Space Station under the following nonlinear initial 


conditions: w-j(O) == 0.1 rad/sec; «2(0) = 0.2 rad/sec.; ^^(O) = 

0.5 rad/sec.; where the desired final state is: = o) 2 ^ = 0; 


- 0.314 rad/sec. = 


It is assumed that three sets of booms 


with a uniformly distributed mass of 1.5 slug/ft can be extended 


along the hub (main part) principal axes from the center of the hub. 
Following the strategy of Ref. 3, and noting that the hub '3' axis 


remains an axis of symmetry during deployment, all booms can be 


extended until = 41.0 secs, at which time the transverse booms 
are stopped and only the extension of the '3' axis booms is continued. 



It is assumed that the principal axis booms have a constant 
extension rate of 4ft/sec. The dynamics of this recovery operation 
is illustrated in Figs. 4.8(a) and (b). 

Suppose that at 50 secs, all operations with the uniformly 
distributed mass booms are halted. At this time the transverse 
booms each have a length of 164ft (49.99m) while - 200 ft 

(60.96m). At t = 50 secs., it should be noted that although 
( 02 /J 2<1 and the nutation angle remains at 17.4°. The composite 

moments of inertia at this time (including .space station plus three 
sets of uniformly distributed mass booms) are: 

= I = 22.91 xlO^ slug-ft2 (30.98 x lO^kg-m^) 

1.3 = 23.82 X 10^ slug-ft^ (32. 204kg-m^) and lo is still greater 
than I. 

At 50 secs, control is initiated using two pairs of offset 
telescoping booms each with end mass of 816kg and an assumed maximum 
length of 100 ft. (30. 48m). For this application both the Q and R 
weighting matrices are selected to be the unit matrix, yielding a 
control law of the form: 

U.,(t) = 0.4996a -1.3249e 
= -0.9998y 

It can be seen from Figs. 4.8 that within 150 secs, of the beginning 
of the recovery operation the nutation angle has been virtually reduced 
to zero and that the x offset boom end mass has a steady state position 
of “22.05 ft (“6.72m ) while the z offset boom exhibits a steady state 
oscillation with an amplitude of 42.7 ft (13.015m). 
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If, on the other hand, we had continued to extend the uniformly 
distributed mass booms along the '3' axis during the same time 
period instead of activating the offset system, even though the 
transverse angular velocity components would have been reduced, 
the nutation angle v/ould have remained essentially constant. 

The effectiveness of the offset system in rapidly reducing both 
transverse rates as well as nutation angle should be noted here. 

I 

Instead of stopping all operations with the uniformly distri- 
buted mass booms at 50 secs., we now consider the case where the 
operations are continued past t = 50 secs, and halted at 75 secs. 

The control is initiated using the z and x offset telescoping booms. 
The time responses of the system as seen from Figs. 4.9 are similar 
in nature to those of Figs. 4,8 except that the time of stopping 
the principal axes booms is changed. At t = 75 secs., it is to be 
noted that I^ is smaller than I and also the transverse angular 
velocities are smaller in magnitude compared to their values at 
t = 50 secs. From the simulation results (Figs. 4.9), it is observed 
that within 200 secs, of the beginning of the recovery operation the 
nutation angle has been reduced to essentially zero. Also, the x 
offset boom end mass has a steady state position of -22.77 ft. which 
is nearly the same value as in Figs. 4.8; but, the z offset boom 
exhibits a steady state oscillation with an amplitude of 16.72 ft. 
which is smaller than the amplitude of 42.7 ft. shown in Figs. 4,8. 
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This is an advantage* since the magnitude of the control effort 


with the z offset boom is reduced when the time of stopping the 
operations of the principal axes booms is increased from 50 secs, 
to 75 secs. It is also to be noted that the control action of 
the offset booms takes place in the linear region of the system 
dynamics (Figs.- 4.9), whereas in Figs. 4.8, the control action of 
the offset booms initially occurs in the quasi-linear region of 
the system dynamics. 
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V. TIME OPTIMAL CONTROL WITH 
SINGLE OFFSET BOOM 


The Investigations presented in Chapter IV considered the 
design of a constant gain regulator for the time invariant system 

X' ^ AX + BU (5.1) 

v/hich minimizes the quadratic performance index. In this chapter, 
the problem of determining the control U (]ul < C*) which forces 
the system (5.1) from the initial state X(n) to zero state in 
minimum time is treated. 

An admissible control U(x), transferring the system state 
from X(0) to X(t^) = 0, is found from the solution of Eq. (5.1) 
given by 

X(t) = e*’' X(0) + r e''*’'"'*’ bu(t) d* ‘ (5.2) 

0 

For X(t.{:) = 0, Eq. (5.2) reduces to 

e'^*^ mi>) del) = -X(0) (5.3) 

0 

As an application of the time optimal theory developed, the 
movable mass system is considered as a two-axis nutation damper 
for the NASA 21 Man Space Station (Fig. 4.2). The standard form 
of the equations of motion, given in Eq. (4.18), is repeated here 
X' = AX + BU 

where 



— 1 

p 

f 


— 1 
o 

-el 


1 — 

o 

It 

X 

3 

,A = 

_d 

1 

o 

,B = 

IJ 


U = n(c"+c) 
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(5.5) 


The quantities d,e,n and the primes have been defined already with 
reference to Eq. (4.18). The solution for U(t), bringing the 
system state to rest in minimum is known to be U(t) = +C*, 
with the number of switches depending upon the initial state of 
the system^ 

Considering the initial states that can be driven to rest in 
a single switch (Fig. 5.1(a)), the control takes the form^^ 

U(t) = K-| for 0 ^ T < Tg; u(t) = for < t < (5.6) 

where 

|K^l=|K2i = C* 

Substitution of Eqs. (5.4) and (5.6) into Eq. (5.3) leads to 

2 

(1 - cosw^Xg) K-j - (cosuqX^ - cosu^Tg) K 2 = a(0)oj^/e (5.7) 

sinco^XgK^ + (sinw^T^ - sinw^x^) ^2 = -3 (0)o)^ (5.8) 

where 

w = /*de 
0 


The expressions for the switching time, x^, and the final time, x^, 
are obtained from the solutions of Eqs. (5.7) and (5.8) as: 




cos 


-1 


(2K2-Ki)K^-(A2+b2) 




■tan~^B^ 

A 


(5.9) 


1 / -1 

= <COS 

f “r ^ 


(2K2-K^)K^+(A2+B2) 


2K2V^A^^ 


■tan“‘‘ ^ 


where 


A'^K-i + ct(0) 0 ) /e 
' 0 

B = 3(0) 0,^ 


(5.10) 
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The control scheme for a single svn’tching for a(0) > 0 
and 3(0) > 0 is shown in Fig. 5.1(b) where 

U(t) = “C* for 0 < T < Tg 

= +C* for Tg < T < . (5.11) 

The system response for an initial condition (a(0)» g(0)) is shov;n 
in Fig. 5.1(c) where the initial state X(0) is driven to rest in 
a single switching at t = t^. The equation of motion of the boom 
end mass is obtained from Eqs. (5.5) and (5.11) as 

^(^"(t) + c(t)) = -C* for 0 < T < T 

" s 

= + C* for Tg < T < 

and hence the time response of the control mass is given by the 
following equation (t < t^): 

C('c) = - -r {u(t) -costu(t)-2u(t”T. ) 
n ^ 

+ 2cos(t-t )u(t-Tc)} (5.12) 

s ^ 

where u(t) is a step input (t>0). 

For the general case, where the initial conditions do not lie 

in the single switching region, piecewise solutions can be used to 

obtain the system response analytically. Also, a time optimal 

control solution can be obtained numerically using the techniques 

17 

of dynamic programming. This approach was employed by Kunciw 
in analyzing the optimal detumbling of the system treated in 
Refs. 4 and 5. 
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VI. CONCLUDING COMMENTS 


As a result of the analysis and numerical results the following 
conclusions regarding the hinged system and the optimal control 
v/ith the telescoping booms can be made: 

1 . Hinged System (Chapter II) 

(a) For the case of a symmetrical spacecraft closed form 
analytical solutions are obtained for two dimensional motion. 

(b) For stability (necessary conditions) certain inequalities 
relating the hinge point offset coordinates to the moment 
of inertia ratio and end masses must be satisfied for the 
three dimensional motion (small amplitude case). 

(c) Hinge damping is always required for the nominal deployment 
of the hinge members. 

(d) Other types of (nutation) damping are requiied to remove 
transverse angular rates effectively since the hinge dampers 
alone do not provide satisfactory nutation decay time constants 
with the system parameters selected here. 

(e) The vertical offset of the hinge points increases the time 
period of the damped transverse rates which may be due to 
the redistribution of moments of inertia for the three 
dimensional motion (small amplitude case), 

(f) The selection of the vertical offset of the hinge points for 
the general case of three dimensional motion for large ampli- 
tudes can not be done easily as there is no criteria which 

I 

assures the stability of the system as in the case of small 
amplitude motion. 
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For the system parameters considered here it is recommended 
that the vertical offset be as small as possible. 

2. Optimal Control (Chapters III, IV and V) 

(a) Three axis linear optimal control of a spinning spacecraft 
system may be achieved by using more than one set of movable 
telescoping booms (with end masses) which are each offset 
from different hub principal axes. 

(b) When only two-axis optimal control is desired only a single 
offset boom is required. The time constants achieved by 
such a system when used for the reduction of nutation 
(wobble) angle are one order of magnitude smaller than those 
previously achieved using non-optimal control logic, 

(c) It is assumed that for the implementation of such a control 
system instantaneous measurements of angular velocity 
components, boom end mass positions, and extension rates 
would be available as well as an on-board computational 
capability. It is hoped to consider the effect of errors 
and delays in measurements on such a system in a subsequent 
study. 

(d) A further application of the offset telescoping boom system 
could be during the terminal phase of a detumbling maneuver 
to quickly remove the residual components of angular velocity 

(e) The time optimal control problem is solved analytically for 
a single offset boom where the initial conditions are such 
that the system can be driven to the equilibrium state with 
a single switching in the bang-bang optimal sequence. 
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VII. RECOMMENDED FUTURE STUDIES 


1 . Further Studies in Optimal Control 

A logical extension of our current application of the linear 
regulator problem for the optimal control of the telescoping system 
would be the examination of the effects of both uncertainties in 
measurements as well as errors in the modeling of the system by 
the linearization technique. The actual implementation of the 
optimal control law, once the gains have been set from a solution 
of the matrix Riccati equation, would involve the sensing of angular 
velocity components by rate gyroscopes and the determination of the 
instantaneous boom lengths, either by a mechanical counter or by 
relating end mass position relative to an optical sensor. In 
other words, physical sensors would actually be involved in the 
real-time determination of state variable components. These 
uncertainties could be included with an application of the estimator 
problem where the differences between the desired state vector 
components aad the actual components due to measurement uncertainties 
would be incorporated within the control logic. The effect of 
errors in the modeling (due to linearization) could be examined by 
simulating the performance of the system using the nonlinear rota- 
tional equations in both linear and nonlinear regions. Finally 
the effect of time delays in the actual measurements (assuming such 
measurements could be performed with zero uncertainties) could be 
analyzed by considering the state estimator-filter problem. 
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The general problem of finding a control law that will result 
in a minimum time recovery could be studied with an application of 
a first-order gradient optimization technique using a digital 
computer. The control variable will involve the boom end mass 
acceleration with respect to body coordinates. Boom position 
can be limited to defined quantities and a penalty function can 
be used to insure a given range of positions. The problem of 
optimal control with a minimum time criterion has been examined 
analytically for the special case of a single offset boom where 
it is assumed that the initial conditions are such that the system 
can be driven to the equilibrium (rest) state with only a single 
switching maneuver in the bang-bang optimal sequence. For this 
system it is possible to obtain an analytical solution for the 
switching and final times in terms of the initial conditions and 
magnitude of the maximum value of the control force, The 
general time optimal control solution for this problem can only 

be obtained numerically using the techniques of dynamic programming. 

17 

This approach has been employed by Kunciw in analyzing the optimal 

detumbling of the system treated in Refs, 4 and 5. 

2. Effect of Gravity-Gradient and Solar Pressure External 
Perturbations 

To date all of the results have been obtained by neglecting any 
and all external perturbations (i.e. in torque - free space). As a 
tumbling spacecraft is recovered by appropriate maneuvers of the 
appendages such effects as those due to gravity-gradient torques and 
solar pressure must be considered. 
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Previous studies have been completed shovnng the effect of 

gravity-gradient torques on the deployment of a gravity-gradient 

stabilized spacecraft which is librating only in the orbit plane 

19 

and where the total libration amplitude is small. Under these 

assumptions a series solution can be used to give an approximation 

19 

of the actual dynamics. Although the present problem can not be 

restricted to small amplitude motions because of the initial tumbling 

or spin, for special cases such as a symmetrical boom deployment 

and motion primarily within the orbit plane it is thought that an 

analytical approach using the techniques of Ref. 19 may be feasible. 

For the general case the effects of gravity-gradient torques could be 

evaluated numerically with particular attention to any readily i den- 

20 

tifiable gravity-gradient induced reasonances. Solar pressure 
torques can be approximated by assuming that for homogeneous boom 
material the booms will bend away from the sun following a circular 
radius of curvature and these torques can be compared with those due 
to the gravity-gradient. 

3. Effect of Flexibility During Boom Deployment 

The problem of including flexibility effects during deployment 

of telescoping appendages has never been fully treated in the open 

literature. A recent examination of a related problem considers 

stability boundaries on the extension of a pair of axial antennas 

21 

whose undeformed state lies along the nominal spin axis. It was 
assumed that the rate of extension was sufficiently small so that 
Coriolis effects due to the rate of change of length could be 
reglected.^^ 
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Ref. 21 determines hov/ far the antennas can be extended before a 
stability boundary is approached but does not simulate the actual 
dynamics during deployment. It is recommended that this problem 
could be considered using two alternative approaches; a) by 
assuming the deployable members to be compound spherical pendulums 
of variable leng.th (according to the control law) where the degree 
of flexibility can be treated by varying effective spring constants 
associated with each of the tvio Euler angles. If deflections are 
assumed to occur primarily in one plane, the spring constant 
associated with the "out-of-plane" Euler angle can be selected 
several orders of magnitude greater than the "in-plane" spring cons- 
stant and b) the booms can be considered as a varying number of 
finite elements depending on how much boom is actually deployed. 

At any given time the combination of rigid body and structural 
dynamics can be examined in a quasi-static manner. The limitations 
of the second method, especially if external continuously acting 
perturbations would be considered later should be noted. However, 
in the absence of external effects it is thought that the tv/o 
approaches could be considered as alternate methods of treating 
the deployment dynamics with first order flexibility effects 
included, and these results then compared with the previous results 
which assume that the booms are competely rigid. 
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COMPUTER PROGRAMS 


1. Hinged System (Chapter II) 

a. Subroutines used: RKSCL, RKGS, SIMQ, RKNXT 

b. Program listing 


IJOB 

1F0BT/A/B/E/P/S fort.ls/l 
IL lSTING 

• c DYNAMICS OF HINGED DEPtOYMENT SYSTEM 

y external HSf?01 r HSP02 

y OI'^EMSTON PARM(5), Yf7), DY(7), ‘a'ORK(F,7), ST7ET7) 

y real 13, T?, 13, L _ 

• COMMON Y 

y COMMON II ,I?,T3,AM,N,M, A,PO,L,CD,RDl,RD2 

y _ COMMON C (E,S) 

y EQUIVALENCE (Y ( 1 . ) , w 1 ) . ( Y ( ?) , W? ) , ( Y ( 3 ) / N 3) , ( Y ( A ) , A 3 ) 

y EQUIVALENCE (Y(S), A?), (Y(61,ALPHA1),(Y(7), ALPHAS) 

y CALL TNOUTCP, S) 

y CALL OPEN (Ir 'SELLAPPAN' ,3,IEP) 

y IFnER.ME. n. STOP UNABLE TO OPEN FILE 

y _ RFAD(?,P1) TMAX,STEP,.T0L , .. 

y ' PEAD(?,Pn SIZE 

y 91 fOPMAT(8F1 0.0) 

y FARM (0 = 0,0 

y PARM(3 )=tmaX 

y ; PARM(3)=8TEP 

y C INITIAL values 

y W I = 0 , 3 

y W2=0,0 

y W3 = a .8? 

y A1 =0 , 0 

• AprO.n 

y ALPHft3=3.F7 

y ALPHA?=I.B7 

y I1=«.R 

y I?=8,S 

y nstn.S 

y AMsO.125 

y RO=i .0 . 

y L = <1.0 

y AsO.O 

y COsO.l 
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RDI = ?.0 
RD?=?,n 
Ms7 
M = S 

WPlTfTC5,9?) TMfiX, STFP, TOL 

lvf<l. A' 'i _ . , _ . ..... - 

KRiTF. rs,9^in. 

WRlTFf'^.Paai ALPHAI, ALPHA? • 

a'Ritf(5^p‘5) t t ,r?, 13 „ . 

WRITF(S,P6) AW.RO.L.A 
W0ITF(5,P7) size 
i'<P!TF(S,9B) 

’p? format/. ' 1 TPAy = ' , FB.?, 1 nx, ' stf.p = ' ,fb. /j , ) ox, ' tol = ' ,fp,6) 

93 format ( '01* Is %F10.B, J ox, 'W?=* ,F1 0,B, 1 OX / '!''^.3= ' ,F1 O.B) 
pai format ( ' OA 1 = ' ,FB. , Rx , ' A?s ' ,FB./4) 
pa? forma T ( ' CAIPHA 1= ' , FA . a,5X , ' ALPHA?= ' , FB.a) 

95 format f 'on = ',FP..a,SX, 'I?z',FP.a,5X, M3=',FA.a) 

96 FORMATC '0AM=',FR,a,5X, 'Rns',FB.a,5X, 'L = ',FR.a,SX, 'As',FA.a) 

97 FORMATC 'OSIZE ' ,7FA.a) 

98 FORMAT ( ' 1 ' , T6, 'T ' , Tl 7, 'K'l * , T30 , 'W? ' , Ta3, 'Vv3 ' ,T5S, ' A 1 ' , ■ 

?T69, • A?' ,TB1 , ' ALPHA 1. ' , TQa, ' ALPHA? ' ,T1 OA, ' IHLF' , T n 5, 'THF TA ' , / ) 

CALL RKSCL (A',SrZF,DY,TOl ,PARM) 

CALL R K r, S ( P A R M , Y , 0 Y , , I H L , H S P 0 1 , H S P 0 ? , A' n R K ) 

WRITF(5,99).THLF 

99 format (' OIHLF s I3I 
CALL EXIT 

FMO 


PROGRAM IS RFLOCATARLF 


TTTL 


main 


FORT ,LS/L 


JFORT/A/B/E/P/S 
i LI ST IMG . 

SURROUTTMF HSROI (T,Y , DY) 

DIMENSION Y(7T, ])Y(71i C(5,5)/ CXf5,5) 
real T t , I?, I3,L 

COMMON wl ,W?,W3, A1 , A?, ALPHA! , ALPHA? 
COMMON I! , I?, 13 , AM,n,M, A,R0,L/CD,RD1 ,RD2 
COMMON C _ „ 

C ' AlsnFP. OF^ALPHAVr A?s0Yr7 OF^ALPHA? 

C CAL. OF COEFF, FOR LHS OF MATRIX EON 

StsSTNC ALPHA n 
S?sSTM ( ALPHA?) 

' Cl sens (ALPHA 1 ) 

c?scos( alpha?) 

A! 1 =?, n* rpn*Rn+A*A+L*L ) 

A1?sRO* (Sl+S?) 

A13rA*(C1 +C?) „ . . 

C(^ , Dsl! +AM*(An*?.0*L*(Al?-Aii)) 
C(i,?)s0.0 
C( 1,3) =0.0 

AiasAM*L#-a + Rn*Sl -A^cn 
c(i ,a)sAia 
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• , (L4Pn*S?-A*C?) 

; C( 1 , S) = At ^ 

; C(?,n = 0.o 

i - . AJ6=?.0*A*A-?,0*A*L4fCl4C?) 

. A17=L4L* (C1*C1 +C?*C?) 

r C(?,2) = I?+A»^4(A16+A17) 

? . A18=A*(Sl-P?)-P0* (Cl-C?) 

f A19=-0.5*L*(?,0*SI*C1-?.0*S?*C?) 

y C(?, •^) = -AK*L*(A1 R4A1 9) 

: C(2,iJ) = o.o . 

: C(?,S)=o.o 

j c ( 3 , n = 0 . 0 

i ... . C (3/ 2K--AM*t4(A) P + Al.9) 

; A20=?.n*Rn*Pn4?.n*RO*L*(Sl+S?)+L*L4r(Sl*Sl+S2*S?) 

; C (3, 3V-T3 + AM^A?0 

t C,(3, ^.0 = 0,0. 

: C(3. 5^0,0 

; C ( ^4 r n A 1 

J c(^»2)=0.0 . . 

; C (a, 3 '1 = 0.0 

J . .... t(.5,A)-M5 ....;. 

; C(5,?)=0.0 

; CC5r3)=0.0 

J c.(5, 443=0.0 

J C(5,S)=AM*t*L 

? DO 1 00 J=1 f M ■ 

; DO 100 1 = 1 

t 100 CX (T » J) = C (T , J), 

t C calculation of PHS of MATRIX' EON 

; 611=12-13 

; BI? = ?,0*A*A-?.0*A>^L*('C1 + C:?) 

; B13=L*l*(CnS{?.n*ALPHAn+C0S(2.0*ALPHA2))-?.0*RO*PO 

• Pli4 = -?.n*P04L*(Sl+S?) 

f Cl V=(H1 1 +AN* (61 ? + Hl 5 + P W4) ) *'«?+w3 

? B15=Rn*(C1*Al4C?*A?)+A*(Sl*Ai4S2*A?) 

? C1? = -?,0*AM^L*fM54WI 

? r3l6=A*CSl-S?)-PO* (C1-C?)-0,5*L>(?.0*S1 *C1-?.0*S?*C2) 

» C 1 .3 = - AN* L * R1 6* ( '/\3 ■* w3-''4?* *•? ■) 

. R17 = A1*M* (RO^C) 4A*Sn-A?4A?*(RO*C?+A*S?) 

• C ia = -AS*l. *617 : : 

t DY ( l)nc 11 +C4 2 + C1 3 + Cl A-RD1*4"1 

; Pta = T3-n 

i P19=L4L*(B1*81+S?*S?)-?.0*(A*A+L*L)+?.0*L*A*(C1+C?) 

; on =(«1 6 + A''’*P 1 O') *i' 

J R?n=A*(Bl -s?)-pn* ( v'1-C?)-0.S* ( 2.0*S1*CI-?.0*S?*C2)*L 

; D1 ? = -AM*L*B?0*l'a 

J B?1=?.0*A*L * (R1*A1 +S2*A2)-L*L*(2, 0*S1*C1 *fl 1 +2* 0*A?*S2*C2) 

• D13 = -AM*r?i _ 

t B22=L4( A1 ^A?)-2,0*A*(ni*Al^C2*A?) 

; P2?0=L*(Cn3{?.uxALPHA1 )*il-COS(2,0*ALPHA?)*A2) 

J .B?21=H2P4R?20 . .... . . . . . .L 

. Oia = -AM*L*P221 *Yi.3 ^ ^ 

• DY(?1=D1 l+DlP + 013 + DVa-RD?*W2 


- 94 - 



B?n = p,0* (■RO*RO + L*L) + ?,0*L*RO* (SI + S?) 
R?5=-L*L*(C1*C1+C?*C?) 

E 1 1 = ( 'A ? + A w * ( R p/i +R p5 ) ) j * i,v p 

P?6 = A*(51-S?V-RO*(Cl-C?)-0.S*L*(?.0*Sl*Cl-2.0.^S? + C?) 
e 1 ?=am*l*p?p^i'‘ 3* *ji 

P27 = ?. n*PO*' (S1*A1 -5^?*A?)-U* (CnS(?. 0*ALPHAn*A1 ) 
8 ?nst*(Al-A?) 4 L*f^nS(?.f)*ALPHA?>)*A?.'> 

E13 = AM*L*(P?7 + B2S)*;^? 

929 = ?,0*RC)*L*(C1*A1+C2*A2) 

830 =L*L*( 2 . 0 *S 1 *C 1 4 ?. 0 *S?*C 2 ) 

E1« = -AM* (b?R4P,30)*w3 
DY(3)=f; 1 t+E1? + En+El<J 
B31=P.0*CJ 40. S4L* 2.0 48 1 *C1 
B3? = A*81-0 .5*L4?. 0*St *r.l 
R33=Rn*Cl +A*S1 . 

B30 = - ( A*C1 +R0*S1-L*C0S (?. 0 4ALPHA1) ) 

E^3yi=cn*i-4L.*Ai 

DY ( y ^ = AM ‘ I, + (R31 *yv3*i^j3 + B32*w?4W2 + R33*wi +‘A't + B3y + iv3*l'*?) -e 30 I 
DY(6) = A 1 

R35=RO*C240, 54t*?. 0*S?*C2 
R38=A*S?^0 ,S*L*2. 04S?*C2 
R37=RO*C?4 A^s? 

f^38 = A4C? + PO*S2-L 4CnS(?.0* ALPHA2) 

R3R1 rCD*L*L*A? 

DY ( = AM4l4 (P3S*W3*W3 + R36 4W2 *w24B37*w14WJ + B3B4'W3*Vj?) -R3Rt 

DY(7)=A? 

call STMCHCX,r)'Y,M,KS) 

IF(KR) 3 ,?f 3 

2 RETlIRK- 

3 WRITE (Pry) 

y EORMAT(//'SIf-Gt'LAR EOUATIOWS') 

RETURN 




PROGRAM IS relocatable 

.TITL HSROl 
IFORT/A/B/E/P/S FORT .US/L 

lUSTING : . : 

• SURROUTINF HSP02CT, Y,DY, IHlF,MnU'V-,P) 

• LOGTCAL RKMXT 

j i Dl'^E.NSJDN Y (7) ,0Y f 7) ,DI)MMY (7) 

dtmfmsion y (7) 

: i PEAl, n ,T2,13,I 

j: COMMON Al ,i-:2jh3, A J , Al. PHAl , A?, ALPHA2 

r : ’ COMMON n ,I?rI3#AM,MrM,A,R0,L,Cn,RDl»PP? 

•i COMMON CISfP) 

ji DEGc‘57.?RS77RS u - ^ - J.- - 

CALL HSRM (T,Y, DUMMY) 

?■ Hi=c( 1 r n*Y(n 

j H2=C(2,?)^Y(2) 

r H3=C(3»3)*Y(3) 

t THeTA=ATAN2(8QRT(Hl4Hl+H2*H2)rH3)*DEG 


- 95 - 



X{1)=Y(U 

X(?)=Yf?^ 

X(5)=Y( 5) 

X ( y ) = Y ( /* ) 

X(5)=Y(5) 

X(6)=YC^)*0FG 
X(7)=Y(7)*nRG 
TP=T+n, 00005 

Ir(.NnT.Ri<r>XT (JHLF)) GO TO ft 

WRTTFf5,1) TP,X’, THLF rTHETA 

1 FORMAT(lX,FO.^f7Pn.7,nO,F9./4) 

WPITF BTMARY(I) T r X ( 3 ) , X C ? ) , X ( 1) , X ( 7 ) , X ( fe ) 

e COMTINUF 

RFTURM- 

ENO 


PROGRAM IS RFLOCAIAFLE 

.TITL HSRO? 


TMAX 

= 10.00 


STFPs 

W1 = 

0,100000 



At = 

0.0000 

A?s 

0,0000 

ALPHA1= 1.6700 


AI.PHA?s 

11 = 

6,5000 

I? = 

P.5000 

AMs 

0,1 P50 

o 

II 

1 .00 00 


SI?F O.JOOO 0,1000 5,0000 


0,05 0 0 Tni, = o , 00 1 ooo 

0,000000 W3= AI,6?0000 

1,6700 

13= '10,5000 

. t"-. » 000 0 A= 0,0000 

5,000 0 5,0 00 0 <14,0000 <11,0 0 0 0 


c, Computer time 

For 10 secs, of dynamic response simulated, 490 secs, of NOVA 
840 computer time were required for a step size of 0.05 secs. 
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d. Subroutine RKNXT 


! JOG 

I FORT/’h/B -'E/P/S FOK r . LO/L 

iLISTIiiO 

; LOuiCFil FUHlTIOH Pi'-f ’.i im>.F < 

; COMHOti ,-R:Kii;-:0- ' 

• LihTh IHOL.il iiE;';T''iP 

] IIiJFF « IHi-F - 

; IF LIP IFF .ER. 0) u-.i 

j IF UDIFi- ^Li’- 0.1 Of 

. McHT “ iF^ 

; GO !U ^'0 

; 35 COr!TiiiLi-u -r-, 

; NtHT = NFXr/.-i5’--i.-IUiF.-.i 

■ 39 CGHTHiOF 

I HOLE iHLF 

: 5 COHTIHjS 

; next = NEXT - 1 

TF LNEXT .or, 0:i GO iO 
; pi<HXT - ,TRiJE. 

; NEXT = 25T<i;IHLF 

GO TO 79 

; 75 CONTI HOE 

; HKHXT = .FhLSE. 

: 79 CONTI MU c. 

; RETURN 

; ENI' 


PROGR.OM TS KF.LOC'^TOPL.E • 


EKiiXT 


ORIGINAL PAGE IS 
OF POOR QUALITY! 
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e. Plotting program 


IJOB 

iFORT/ft/R/f-:/»/S FOPT.IS/L 
i LI STING 

j PARAMFTFR NPTS = ?01. , NV=S 

• DIMENSION T(NPTS), VCNPIS^NVI 

? c 

? CALL TNnuT(?,F) 

j CALL SPFFL("!^FLT% lER) 

; IF C'UR .NF, }) vSTOP -- ON SPOOLING 

? WRITF(S,00) 

? CALL OPENdr 'SELLAPPAM' ,1.IEP) 

? IFCTER.NE.n ■STOP unable TO OPEN FILE 

; READ binary ( 1 ) ( T ( T ) , ( V C I , J) » vJ = 1 , N V ) , I 

• DO A M=1,NV 

t REAOf.?f-n WIDTH, height, T'YA y^VMJN^V^^Ay 

• WRITF(5,9?) WIDTH, HEIGHT, TMAy, VMIN, 

• ' CALL PSI7F. ('WIDTH, HEIGHT) 

: CALL (^pny 

; CALL PAXES 

• CALL PLOPCP.O,T,Tiw-,Ay,\/Y'lN,vn ,M),VMAy, 

• CALL PLOTCl .0+HEIGHT,0.0,-5) 

f 8 CONTINUE 

; CALI EXIT 

. 9 FnRMATCBFIO.d) 

• 9D FORMATCIM 

• 92 FORMAT (2X, /, IX, 1P5FIE>.5) 

• END 


PROGRAIM IS relocatable 

•TITL .MAIN 
JRLDR/H TMP/S 001 DPOrPLOT.LR FORT.LH 
lEXEC 


=1 ,NPTS) 
VMAX 

NPTS) 
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2 . Optimal Control (Chapters III and IV) 

a. Subroutine used: HINV, RKSCL, RKGS, SIMQ, RKNXT 

b. Program listing 

PROGRAM I 


I JOB 

iFORT/A/B/E/P/S FORT.LS/L 
lUSTlNG 

; C PROGRAM I 

•' C PROGRAM TO FIND A AND B FROM G-H-F MATRICES 

t D I MENS I ON A c 3, 3 1 , B ( 5 , 2 1 , G C 5 . 3 ) » F ( 3 / 5 ) » H ( 3 , 3 ) , F ( 3 f ? ) r 

• *LL(3),mm(3) 

; REAL II, 12, 13fMK,Mi ,M2,LM 

• CALL INOUT(?,S) 

: 106 FORMAT (5X , 13H THt A MATRIX/) 

y 107 format (5X , 1 3H I HE B MATRIX/) 

t lOB FORMAT(6(-1PF?O.B) ) 

j N=3 

? L = 2 

; RE AO (2, 101) 11/ 12, 13, MR 

; WR1TE(5, 1 01) 1 1 , 12,T3,MH 

; 101 format (XlF 2 0,2) 

; REAOr?, 102) Ml ,M2,LM 

? WRITECS, 1 02)Mt ,M2 ,um 

; "102 FORMAT (aF20. 2) 

? RtAD(2, 103) AA,HB,r.C 

T WHITFC5, 103) AA,BB,CC 

j 103 format (iiF20. 2) 

: Cl =AA/LM 

j C2 = BR/L'M 

V* •' '"C3 = CC/LM ■■ ■■': . , 

t U2=M2/CMb+M2) 

• U = Mlyt-(MB + M?) /(MH + Ml +M?) 

f ULM = U*l-M*LM 

; AI1=T 1 /ULM 

t A12=I?/ULM 

; AI 3sT3/OLM 

; C CAL. OF G MATRIX 

? G(l,l)=AIl+C2*C2+C3*C3-2.0*U2*C2^C3 

) ~ ' G(l ,2)=-Cl*C2 + U2*Cl*C3 

t G(l,5)=0.0 

t G(2, 1)=-C1*C2+U2*C1*C5 

ORIGINAL PAGE IS 
OF POOR QUALITY 
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t 


G(2 

, ?) = AT?+Ci*Cl 

• 

r 


G(2 

,3) =0.0 



G(3 

, n=o.o 

• 

f 


G(3 

,2)=0.0 

• 

9 


6(3 

, 3) = AI3 + C) *CH'C2*C2 + C3*C3-2.0*U2*C?*C3 

9 

9 

t 

cal 

,df h matrix 

• 

r 

\ 

H(1 

, 1 )=-CUC?-t'U2*Cl*C5 

• 

9 


H(1 

, 2)=-( AI 3-AI2 + C2*C2 + C J*C3) +2.0*U?*C?*C3 

• 

9 


H(1 

,3)=0.0 

! 


H(2 

, 1 )=-(Aii-Ai3-ci*cn 

9 

9 


H(2 

, 2)=CI*C?-U2*C1*C3 

t 


H(2 

,3)=0.0 

9 

9 


H(3 

, 1 )=0.0 

• 

f 

- 

H(3 

,2)=0.0 

9 

9 

! 

H(3 

,3) = 0.0 

9 

9 

C 

cal 

. of f matrix 

9 

9 

i 

F(1 

, 1 ) =-C2+U?^C3 

9 

9 


F(1 

, 2 ) = 0 . 0 

9 

9 


F(2 

fl)=Cl 

• 

9 


P (2 

,2) =0.0 

9 


F(3 

, n = 0 .0 . 

9 

9 


F(3 

,2)=C3-U2*C2 

9 

9 

c 

inverse of g=e 

9 

f-. 


DO 

13 1 = 1, N : 

• 

9 


DO 

13 J=1,N 


; 13 £(1, J)=r,(I , J) 

• CALL MINV(E ,N,D,LL»MM) 

V"C“’ CALrOF A = E*H 

; WHITE( 5 » 10 ^O 

. DO IS 1 = UN 

• 00 15 J= UN 

• A(UJ) = 0.0 

DO IS K = UN 

■; is' AUf J)sA(U J)+E(I »K)*H(K, J) 

j DO 16 1 = UN 

? 16 i^jHITEC5,108) (A(UJ),J=UM) 

• C CAL. OF B=E*F 

t WRITE( 5 » 107 ) 

DO 17 1 = UN 
j “ DO 17 J = UL 

; 6 (UJ)= 0,0 

; DO 17 K= UN 

17 5 ( 1 , J)=R(UJ)+ECUK)*F(UJ) 
t DO I B I = U N 


18 v^RITF (5, 1 08) (B(U J) »0 = UL) 

call exit 
end 



program is relocatable 


,TITL .MAIN 

IRLOR/M TMP/S’ OOi DPOtSSP.LB FORT. LB 
lEXEC 


-TOO- 


10500000, 

00 

1250000'^* 

00 

15000000, 0< 137000.00 

1800. 

0 0 

1800. 

00 . 

100. Of 

65. 

00 

0. 

00 

65.0 0 

THE A 





2.161900006 

-3 

-5.55732^006 

-1 

O.OOOOOOOOE 0 

6.001A19006 


-2. 161896006 

-3 

O.OOOOOOOOE 0 

0.000000006 
THE 8 MATRIX 

0 

0.000000006 

0 

0.000000006 0 

8,277688006 

"3 

0,000000006 

0 

■ '■ 

5.7/il83700E 

-1 

0.000000006 

0 


0,000000006 

0 

3.808111006 

- 1 



PROGRAM II 


I JOB 

iFORT>A/B/E/P/S FORT.L.S/L 
aiSTING 

; C PROGRAM II 

? C' KICCATI FQtJATIQN' PROGRAM CRTCATI) 

t OIMFNSlOf^ A ( 1 0» i 0) rFU 1 0, in ) ,0 (1 0, I 0) , R (1 0, 1 0) 

t * K( 1 n, 1 0) ,0(10, 10) (10) ,F n 0, 1 0) ,E( 1 0, 1 0) r 

■ ? ■ ; V G( 1 0, 10) ,H ( I 0, 10) , S( 10, 1 0) 

?■ REAL K 

j INTEGER OPTION, BLANK 

: ""''COMMON/KALKAN/ICC, IFF, blank 
; DATA ICC, IFF,BLANK/*C ','F V 

; CALL TN0UT(2,5) 

"r 10 0 0 FORMAT ( IHl ,SX, 37H0PT I MAL CONTROL/) 

i 1001 F0RMATC2I?) 

r 1003 format ( IHO, 5X, 13H the a matrix/) 

100^ format ( 1 HO, 5X , 1 3H THE B MATRIX/) 

? 1006 FORMAT(8F10.3) 

; 1 007 format (Al ,9X, RF10.3, 13) 

r008 '‘F0RMAT ( 1H0,^B(1H*) ) 

; 1010 format ( 1H0,SX,P2H*** CONTROL OPTION ***/) 

r 1011 F0RMAT(6( IPEPO.B) ) 

V lOia FORMAT(IHO,BX,13H the R matrix/) 
r 1013 FQRMATt 1H0,SX, 13H THE Q MATRIX/) 

; 101 « FORMAT (1H'0,5X,19H INITIAL CONDITIONS/) 

I"' iOlS format C1H0,5X, 8H TIME = , 1 PERO . 8/6 X , 5HG A INS ) 

; 1016 formatciho,5x,21hsteady state solution// 

; * 6X,6H GAINS/)' 

• 100 RtAD(R,100l) N,M 

i write (5, 1000) 

t _ WRITE C5, lOOn NfM 

"r"" I T E C 5 , 1 0 0 6 ) ■ ^ 

? WRITECB, 1 003) 

j DQ no I = 1,N 

• R£AD(?,1006) (ACi, J) , J=1 ,M) 

f no. WRITE (5, 1 01 1 ) (A(I , J) , J=1 ,N) 
r WRITEC5, lOOA) 

I DO 1?0 m,M ' 

t REAfXpn 006) (H( J, n f J = 1 f 

? 120 WRITEC5, lOin (B(J,I),J= 1 ,N) 

t ISO READ(?, in07,END=R99) OPIION , T 1 , T ? , NPT 
; WRITE (B, 1 006) 

; IF(0PTI(1M.F0,PLAMK) GO TO 100 

J IFUiPnON.EO.ICC) GO TO 500 

? 300 WRITE (5, 1 010) 

f NRsM 

..... 

t DO 330 1=1 ,N 

• DO 310 J=1,M 

■"Y"- ‘510 Ea,j)=B(i,j) 

; DO 330- J = 1 ,N 

) 330 F(I, J)=A(J,1) 


- 102 - 



GO TO a 0 o' 

400^*'H1TE(5#1012) 

DO a 1 0 1 = 1 » nr 

READ(?fl 006) (R(I » J) / J = l 

aio WHITE(5,10in (R(IrJ)»J=l»NR) 
ARITF (5, 1 0T3) 

DO U?0 I = l,NQ 

READ(2,1006) (Q( 1# J) f J=1 ,N0) 
^i20 WRITEfSMOH) J=l,MO) 

DO ajo I = WNR 
DO a30 J=1,NR 
aso H(i, j)=R(i,j) 

”” DO 4410 1 = 1 , NR 
DO 44440 J=1/N , 

R(I»J)=0.0 
DO aoo II = 1,I\!R 

4400 R(j , j)=Rri,j)+H(i, in*E(j, in 

IK(OPTIOM.EO.TCC) GO TO 500 
500 'IF(NPT.GT.O) GO TO 530 
DO 5?0 I=1,W 
DO 510 J=1,N 
5li0 G(nj)=0.0 
520 G(I,I)=1.0 
_ tPS=0. 1 
go' to 570 

530 WRITEC5, 1 0 14)) ' 

DO 54J0 I = 1,N 

READ( 2,1006) (Gri, J) , J=1 fN) 

544 0 WHITE (5r 1011 ) {G(I, J), vI=l,N) 
WRITE C5, 1 0f4B) 

TlMe=ABS(TP-Tl) 

PTS=200. 0*TIME 

XX = NPT 

XXsPTS/XX 

1 D= X X 
DI=ID 

YY1=ABS(XX-DI ) 

IF( YY,GT,0,05) ID=ID+! 

ITsPTS 

EPSs'o . 005 

IF(OPTION.EQ.ICC) TIME=T2 
WHITE (5, 1 015) TIME 
DO 560 1=1, NR 
DO 550 0=1, N 
Kn,J)=0.0 
DO 550 1 1 = 1 ,N 

550 K(I, J)=K(I ,J)+R(I,II)*G(II, J) 
560 WRITE (5, 1011) (KCI, J) , J=1 ,N) 
“570 LC=0 

ICHsIO 

575 DO 580 1 = 1 ,NP 
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r i ; 00 J=i,N 

t\ K(I,J) = 0.0 

r DO spo u = i 

r sao K(I ,JY=K (I , J)+R( I , II J) 

J DO 590 1=1 ,N 

r DO 590 J=)',N 

r H(i,j)=o.o 

t DO 590 T 1=1 ,NR 

f 5^0 

'ii " i 00 ^10 1=1, N 

% no 610 0 = 1 , N 

? D(I,J)=Q(I,J) 

j DO 600 11=1 ,N 

i 600 D(i,j)=i){i,Jv+F(i,ii)*G(ii,j)+Gn,.n)*F(j,in-G(i,ii)*HnT,J) 

t 610 S(I, Jl=Gn,J)+D(I,J)*:EPS 

IF(NPT.LF.O) go to 640 

; LC=LC+1 

: ' IFac„LT.ICHl GO TO 625 

f ' 1CH=TCH+JD 

f ALC=LC 

; : TsALC*EPS 

! ' IF(OPTIOO.FQ. ICC) TlME=T2-T 

? WRITE (5, 1015) TlMfc 

J : no 620 1 = 1 ,WR 

r ' '62 0 ■ WRITE (5, 101 1 ) (K(I,J) , J=1 ,N) 

? IF (LC. OF. IT) GO TO 150 

f 625 DO 630 1=1 ,fvl ' 

f . DO 630 J=1 ,N 

? 630 G(I,J)=S(I,J) 

J i GO TO 575 

J 640 SU«=0, 0 

t DO 650 1=1 ,N 

; ; i DO 650 J=l ,N 

f i SUMsSUM+ABSCDCl, J)) 

J 650 6(1, J)=S( I , J) 

? IFtSU'M.GT.O.Ol ) GO TO 575 

WRITE (5, 1016) 

} DU 660 I=l ,NR 

J 660 WRITE( 5 , 10 in (K ( I , J ) , J= 1 , N ) 

1 i i GO TO 150 
; i 99:9 CONTIMUF 

f CALL EXIT 

; ' ” ENO 


PROGRAM IS RELOCATARLt 

,TITL .main 
IHL0R/_M TMP/S OOl FORT. LB 
'I EXEC" 


ORIGINAL 
OF POOR QUALITY 
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OPTIMAL CONTROL 

i'?r 


THE A MATRIX 

2. 1 6190000 E ' - 7 > 
6.00191900E -1 
p.OOOOOOOOE 0 

THE R MATRIX 

e.27769?00E -3 
O.OOOOOOOOE 0 


*?.*,. CONTROL OPTION 


„the; R matrix 

1 . OOOOO’OOOE, 0 
O.OOOOOOOOE 0 

THE Q MATR'IX 


9.B6000100E 0 
O.OOOOOOOOE 0 
O .OOO OOOOOE _ 0 


STEADY state SOLUTl 


GAINS 


-2.21919500E 

' "'O.OOOOOOOOE 
iEOF 


0 

0 



105 


-5.55732^i00E -1 'O.OOOOOOOOE 0 

-2.16190000E -3 O.OOOOOOOOE 0 

O.OOOOOOOOE 0 O.OOOOOOOOE 0 


5.7A183700E -1 O.OOOOOOOOE 0 

O.OOOOOOOOE 0 3.8A0in00E -1 




O.OOOOOOOOE 0 
l.OOOOOOOOE 0 


O.OOOOOOOOE 0 
9.86000100E 0 
O.OOOOOOOOE 0 


O.OOOOOOOOE 0 
O.OOOOOOOOE 0 
9,860001 OOE 0 


3,e0596200£ 0 O.OOOOOOOOE 0 

O.OOOOOOOOE 0 3.19006300E 0 


PROGRAM III 


IJOB 

iPORT/A/B/e/P/S FORT.LS/L 
IL I STING 

C PROGRAM in 

C ' optimal control « TIME RESPONSES 

external BSOCWRSQC? 

DIMENSION PARMC5), Y(7),OYC7),WORK(R,7)fSI2E(V) 

REAL n , I?, ,M2,L''‘- 

; HEAL KCM f KC12,KC13^KC21 /KC2?,KC?3 
COMMON Y 

! : common n, I2r I3,MB,Mt,M2,LM 

' ^ “ COMMON N,M 

' i COMMON AA,BBfCC,00 

j COMMON KCU»Ka?,KClS,KC21 /KC?2,KC23 

t ' ; ' r" COMMON CC3f 3) 

j equivalence (Y (U , WIN) , ( Y(2) , W?N) , (Y(3) »rt'3N) , (Y (U) ,DZN) 

; equivalence (Y(5),DXN),(YC6),ZN), (Y(7).XN) 

Ji : ^ ' CALL INOUT(?,E) 

?: call open (1, 'SELLAPPAN* ,5fIER) 

; IFdER.NE.n STOP UNABLE TO OPEN FILE 

; 'READ (2, 91) TMAX,STEPdOL 

f READ(2r9l) SIZE 

j: 91 FORMA T(8F 10.0) 

)' "‘PARMn )s0.0 

j PARM(?)=TMAX 

t i PARMC3)=SieP 

t ■initial values 

• READ(2,9l 1 ). I 1 , I?r I3,MR 

: READ (2,911) M1,M2,LM 

' "■■REA0(2,91 n AA,BB,CC,DD 
t 911 FORMAT (UK 20. ^) 

?| WlNsO, 0391/0. 31« 

> ' W2N = 0,0 

j| W3N=(0,35-0.31 A)/0,3ia 

yi _ DZN = 0,0 

> 'DXNsO, 0 

i ZN=0,0 

S REA0(2,Q15) Kcn ,KC12,KC13 OF^PO^ 

j READ (2,91 3) KC21 ,kC22,KC23 QUALITY 

• 913 FORMAT(RFIO.n) 

i N = 7 

j ; '-'=3 

y ' WRITE(5,9?) TMAX,SIFP,T0L 

y : ARI1E(5,95) 11,12,13 

y WRITE(S,P6) M1,M?,LM 

y _ i wRITE(S,97) AA,B8,CC,DD _ „ 

F" write (S, 98) SIZE ■ ■ 

y WRITE(B,9R1) KC 11 , KC I 2 , KC 1 3 

y wRlTE(S,98i) kC21,KC?2,kc?3 
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; 92 FOWMATf 'nMAx=»,p8.a,10>:, 'vSTEP=',FB.a, lOX, 'TnU = *rFP.6) 

; 9S FOPMAT('OI1 = SF1S.2,2X, 'TP='rFlS.2,?X, M3sSF\5.?) 

t 96 FOWMAT('OM1=',F15.?,2X, 'M2r',Fl5.2,2X, »LF=',F1S.2) 

; 97 format ( *0AA = ',Fp.q,5X,'RR=',PB.a,5X,*CC=%Fe.a,SX, '00=' ,F*8.a) 

f 98 format ( 'OSIZE' , 7F6.a) 

• 981 format ( 8F1 0. a) 

• 99 FORMAT ( ' 1 S T6, M , T1 7f 'M S T30^ 'W?', T^3» ' ,T88, 'DZ » , 

i 2T69, *DX',TP1 , »Z»,T9a, 'X',T10R, 'THETA',T115, 'IHLFS/) 

; . CALL RKSCLCf'i, SIZE, DY, TOL, FARM) 

J CALL RKGS(PAPM,Y,DY,N,XHLFr8SOCl,BSOC?,i'JOPK) 

; . WRirF(SrlOO) IHLF 

j ioO FORMATC '.OIHLF=S 13) 

J CALL EXIT 

; END 




PROGRAM IS RELOCATAF^LE 

.titl .main 
IFOR r/A/B/E/P/sS FORT.LS/L 
ILISTING 

: SUBROUTINE BSOCl CT, Y,DY) 

I dimension Y(7) ,0Y (7) ,C(3, 3) 

REAL II, J?, I.3,MB,M1 

REAL KCn ,KC12,KC13,KC?1 ,KC2?,KC?3 
common VM 1 N , W?N , W3N , DZN , DXN , ZN , XN 
COMMON I 1 , I 2 , I 3 , MB , M 1 , M? , LM 
COMMON N^M 

common aa,hb,cc,od 

COMMON KC11,KC12,KC13,KC21,KC22,KC23 
COMMON C 

CAL.OF.COEFF.FOR LF-iS OF MATRIX EQN 
U=M1 * (MR+M2) / (MB+M1 +M2) ^ 

■■’ULM = LULM»kLM ' - "■ 

NON 0 I MENS I ON A L 1 Z A T I ON 
i Alt^n/ULM 
AI2=I2/ULM 
? : AI3S13/ULM 


C1=AA/LM 

C2SBB/LM" - 

C3=CC/LM 
C9=DD/LM 

Ul=‘M * (MR + Mp ) / ( (MR + Mi +M2) *in 
U2=M2* (MR+Ml)/ ( (MB+Ml +M2 ) *U) 

u i = *Mp/ ( + N J +wp) *D) 

C ( I , n=A II +U1 * (C2*C2 + ZN*ZN) +U2*(C5*C3 + C9*Ca) 
1 : ; *-2,0*U3*(C2*C3+C9*ZN) 

;.J_: : C(l,2)=-U1 *(CI*C2)-U2*C3*XN+U3*(C1*C3 + C2*XN) 

? C(1 ,3)=-Ul*(Cl*ZN)-U2*C^*XN + U3*(Cl*Ca + ZN)kXN) 

J CC2, 1 )=*U1*(C) *C2)-U2*C3^XN+U3*(C1*C3+C2*XN) 

? C (2, ?) = AI2 + IJ1*(CI *-Cl+ZN*ZN)+U2*(CU*Ca + XN*XN) 

; . *-2.0*U3*(CI*XN+C9*ZN) 
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C(2,3)=«U1*(C?*7N) -U 2*C <i + U3*(C2*Ca + C3*ZM) 

C(3» 1 )="U1 *{C1 *Z(\!)-U2^Cil*XIM + U3*(Cl*C/i + 7.N*XN) 

C(3, 2)=”UJ *CC?*ZN)-U2*C3*C«+U3* (C2*C/4 + C3*7N) 

C(3#3)=AI3 + U1 *(C5 *Cl+C2*C?)+U2*CC3*C3 + XN*XfvJD 
*-?.0*U3*(C2*C3+C1*XN) 

cal. of RMS OF J^-ATRIX FON 

AZN = - (KCli 1 M+KCi ?*In2N + KC15*W3N)-ZIv 

IF (M2, to, 0.0) GO TO 7b 
F1=2.0*U3*C1/(U2*C3-U3*C2) 

GO TO 76 

75 F 1=0,0 

76 AXN = "(KC21 *W1N + KC22*'<^2N + KC23*w3N)-F1 *DXN 
SM = 1 . 0 + ■•■\'3^ 

A11 = (AT3*’A12)*w?N*SN 

V A12l=-C1^ZN^W1M*W?NJ + (C2*C?-ZO*7N)*w2N*SM + Cl*C2*WlN*SI'i 
A122 = ?.0*ZN*DZO*‘^'1M + C2*ZM* (SN*S0-W2N*W2N)+C2*A7fJ 
' A 1 23 = “C y * XR* 1 0* w20 + ( C 3*C3"C 0 *C 0 ) *W?N*SR + C3* XM* w 1 N*S(\! 

!' A120 = -2, 0*C3*DXL>i'J2N-2.0*CO*DXM*^SiM + C3*C^*(8G*SN-W2N*W2N) 

; A1? = U1*( A)2l+A122)+LI2*(A1?3 + Al2a) 

^_A15 = " (C 1 *Ca4‘ZNycXN) *w lN*W2N+2, On (C2*C3-C^J*ZN ) nwpNnSN 
: *A1 6= ( C 1 *C3 + C2*XN ) nwlN*SM + 2 . 0*CO*nzOniMN-2, 0*C2 *DXN*‘a 2M 
: AJ 7 = -2. onDXM*7 0*SN + (C2*C0 + C3*ZN) *(SN*v8N-i'.‘2N*W2l'j) + C3*AZN 
A13=U3*(A55+A1 6+A17) 
bY(l)=-(All+Al2-A13) 

B1 1 = (AI l-A,I3)nSM*l‘.lN 

,, B121 =-C2nZNni\tfvJ*i\2N + C 1 nC2*W2N*Sf'J+(Cl *C1~70*ZH) *SN'*W1N 
B122 = -2, OnZOnOZNnwpN + C 1*ZM* (SNnSN-wi N*w 1N‘) + C 1 * AZO 
B123 = «'C3nCa*WlK'nv'/20 + C3nXN*W2N*SM + (-C^*CA + y Mnx'iv) *SOnwl M 
: B12^ = '-?:,n*XM**'XN*W2M + C0*XO*(r.NnSK!-WlN*rtlM)-C0nAXN 
Bt 2 = U1 * (R121 +R1 22) 4-02* (R) ?3 + ftl?0) 

B15 = -(C?*C^44C3*ZN)*w1M*i'v 2W + (Cl *C3-!-C2*XM)*l\'?N*Sf-* 
i_Rl6=?,0*(Cl *XN-C^*ZN) *SN*'AI10-2 , 0* (C1*DXN+C«*DZM) *W2N 
' B17 = (C1*C/4 + XN*ZN)*(,SN*SM-w1N*w1N)+X(m*AZN-ZO*AX!\i 
B 15 = U3* (R154-B1 6 + Bl 7) 

OY (2)="(B1 1-B124-B13) 

. C11 = (AI2-AI 1 )*WliM*w2N 

' C12l = (C2*C2-C 1*C 1 ) *wif\i*K'2N-Cl *ZM*W2M*SN4'C?*ZN*S0*W1M 
[ C122 = 2, 0*C l*nZN*v‘,'li'J + 2. 0*C2*D7N*1'^?N4'C1 *C2*( W1 M*Wl N-W?N*W2M) 

; C123 = (C3*C3-XN*X.M) *i\' 1 M*W2M«C0*XN*W2AJ*SN + C3*C A*SN*W1 N 

' C12<4 = -2,0*X^J*nXRi*SO + C5*X'M* (W1 N*’a'1NJ»Im 2M*W2N) + C3*AXN 
C12 = U1 *(C121 +C12?)+U2*(Cl?34-C12a) 

C 1 5 = 2 . 0 * ( C 2 * C 3 "C 1 * X Iv ) *l'J 1 N * W?N - ( C 1 *C a 4- X R * 7 N ) *1‘‘2 M * SN 
C 1 6 = ( C 2 * c di + c 3 * Z M * S .M * A 1 M 4- 2 . 0 * X M * D 7 N- * *V 1 N 4 ? , 0 * C 3 * I '> X N * 2 N 
C17 = -2,0*Cl*nXM*SM4(Cl *C34C2*XN) * (I'M N + ai1N-w 2M*1-'2N) 4C2* AXN) 
Cl3 = i.)3* (C1R4C1 64C1 7) 
l)Y(3)=-(Cl 1-C124C1 3) 

DY (dO sAZN 
"■'DY (5) = AXM 
DY(6)=DZN 
_DYCM=OXO 

CALL SIMQ(C,0Y»M,KS) 

IF(KS) 3,2,3 




D5' PdGR 


2 WfcTURM • 

3 ivKITt(5,a) 

tx FORMATC// ' SINGULAR EOUATIONSM 
RETURN 

END ■ ' ■ 


PROGRAM I's”RELOCATABLE 

.TITL BSOCl 
IFORT/A/R/t/P/S FORT.LS/L 
'ILISTING 

y SUaROlJTINt PSOCS ( T , Y , 0 Y , I HL F , NDUM , P ) 

t i LOGICAL RKNXt 
y" DIMENSION Y ( 7 ) » 0 Y ( 7 ) f OIJMM Y ( 7 ) 

y dimension yC7)- 

y REAL 

f REAL KCl r,KC12.KCl 3,KC21 #KC22,KC23 

common vMN,w2N,N3N,DZN,DXN,ZM, xn 

COMMON II, I?,I3fMR,MlrM2,LM 

common N,M 
common aa,bb,cc,od 

COMMON Kcn ,KC)2,KC13rKC2l ,KC22,KC23 
DEG=57.?R577P5 
CALL RSOCUTrY, DUMMY) 

WS = 0.31^i 
xcnsYcu*ws 
X(2)=Y(2)*i‘'S 
X(3) = ( 1 .0 + Y(3) )*'^'S 
X(a)=Y(y)*LM*ws 

X(5)=Y(B)*LM*v.‘S 

X(6)=Y(6) *LM _ , __ 

XI7)sY(7)*LM’ 'T r 

Hi=ii*x(n " 

H2=I2*X(2) 

H3=15*X(3) 

THETA = ATAN?(SQPT (H1*H1-«-H2*H2) ,H3)*0EG 
lAsT/wS 

T P = r A + 0.00 005 ' ■ ■ 

IF(.NOT.Rknx7(IHLF)) go to 8 
y! : WRITE(S,1) IP, X, THETA, IHLF 

y. : 1 FORMATClX,FR.a,7F13,7,F9.y,no) 

y kMRITE BINARY(l) T A , T HE T A , X ( 7 ) , X ( 6 ) , X ( 3 ) , X ( 2 ) , X ( 1) 

y 8 CONTINUE 

y RETURN 

y end 


PROGRAM IS RELOCATABLE 


TITL fiSOC? 


TMAX= 62.fi000 


STEPS 

0.3100 


T0L=0. 001000 

11 = ' ib^>ooooo 

.00 

12= 12500000.00 

13 = 

15000000.00 

Mir. 1800 

o 

o 

M? = 

1800.00 

LMn 

100.00 

AA= 65.0000 

HR = 

0.0000 

CC= 65, 

0000 

DDs 0,0000 

SIZE 0. 1000 0 

-2.219? 3. 

0,0000 0 , 

. 1000 
8060 
00 00 

0. 1000 
0,0000 
3,1001 

0,0000 0 

, 0000 

2.0000 2,0000 


c. Computer time (Optimal Control) 

PROGRAM I: 109 secs, (matrix inversion) 

PROGRAM II; 120 secs, (to reach steady state) 

PROGRAM III: 382 secs, (to simulate 200 secs, real 

time response with stpe size= 1 sec.) 
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